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Introduction 

The  scientific  objective  of  this  research  was  to  develop  and  characterize  an  MRI-coupled  fluorescence 
molecular  tomography  (FMT)  system  and  associated  reconstruction  algorithms  to  image  exogenous  fluorescent 
agents  in  human  breast.  Imaging  the  spatial  distribution  of  unspecific  and/or  targeted  fluorescent  probes  may 
help  determine  whether  a  suspect  region  in  an  MRI  image  is  indeed  malignant,  thus  reducing  the  high  number 
of  false  positives  and  associated  biopsies,  and  aiding  in  treatment  planning.  Fluorescence  imaging  at  depth  in 
tissue  is  challenging  for  a  number  of  reasons.  The  nature  of  non-invasive  diffuse  imaging  requires  detection  at 
the  tissue  surface,  making  the  data  set  fairly  sparse  compared  with  the  volume  of  interest.  This  results  in 
relatively  poor  resolution,  depth  dependent  sensitivity,  and  inaccurate  quantification,  especially  in  complex 
tissue  volumes  such  as  human  breast.  Secondly,  low  levels  of  fluorescence  emission  are  often  swamped  by 
excitation  cross-talk  from  the  exciting  laser  source.  The  FMT  system  here  seeks  to  address  both  of  these  issues. 
Image  resolution  and  quantification  performance  is  addressed  by  incorporating  the  system  into  a  clinical  MRI 
for  dual-modality  image  acquisition  and  using  the  MR  image  as  a  tissue  structure  template  to  guide  the  FMT 
image  formation.  Excitation  contamination  is  addressed  by  using  spectrometer-based  detection  systems  which 
provide  a  unique  opportunity  to  de-couple  fluorescence  and  excitation  signals  using  a  spectral  fitting  routine. 
This  report  summarizes  the  major  accomplishments  completed  during  this  fellowship. 

Body 

The  work  completed  for  this  traineeship  accomplished  the  proposed  aims. 

Aim  1:  It  was  shown  that  the  system  is  appropriately  sensitive  for  FMT  imaging  during  standard  breast-MR 
protocols. 

Aim  2:  All  results  from  numerical  and  phantom  data  showed  dramatic  improvements  in  imaging  performance 
with  MRI  prior  information  incorporated  into  the  imaging  algorithm.  In  most  cases,  the  MR  information  was 
required  for  accurate  imaging  and  the  image  reconstruction  problem  was  intractable  without  the  MRI  images. 
Aim  3:  Numerical  and  phantom  data  was  used  demonstrate  FMT  imaging  in  heterogeneous  domains.  It  was 
shown  that  the  MRI  information  is  required  to  produce  meaningful  images  in  these  complex,  though  realistic 
volumes. 

Aim  4:  This  aim  was  revised  when  Pharmacyclics  halted  manufacturing  of  LuTex.  The  fluorescent  drug  used 
in  these  studies  was  transitioned  to  Indocyanine  green  (ICG),  a  non-specific  fluorophore  approved  for  use  in 
humans.  It  was  shown  that  ICG  can  be  quantified  down  to  at  least  1  nM,  much  lower  than  the  50  pM 
concentrations  required  for  most  gadolinium-based  MRI  contrast  agents. 

Aim  5:  It  was  shown  that  the  signal  producing  the  most  sensitive  measurement  depends  on  the  concentration  of 
the  drug.  As  drug  concentration  increases,  the  signal  providing  the  most  sensitive  measurement  transitions 
from  fluorescence  to  absorption. 

Imaging  System 

One  of  the  primary  accomplishments  of  this  project  was  the  design  and  development  of  a  parallel  spectrometer- 
based  tomographic  imaging  system  which  couples  into  a  Philips  3T  MRI  magnet.  The  system  was  developed  to 
acquire  spectrally  resolved  transmission,  fluorescence  emission,  and  bioluminescence  spectra.  Photographs  and 
a  diagram  of  the  system  are  presented  in  Fig.  1.  The  spectroscopy  system  is  composed  of  sixteen  Insight 
spectrometers  (Acton  Research)  with  Pixis  cooled  CCD  cameras  which  provide  highly  sensitive  spectrally 
resolved  detection  of  the  collected  light.  Each  spectrometer  is  coupled  to  the  patient  interface  via  a  1 3  meter 
bifurcated  silica  fiber  (a  total  16  fibers)  custom  designed  for  this  system.  These  custom  designed  bifurcated 
fibers  circumscribe  the  tissue  volume  being  imaged  and  are  used  for  both  source  delivery  and  detection.  One 
end  of  each  bifurcated  fiber  is  coupled  into  a  spectrometer  through  a  custom  designed  and  manufactured  input 
optical  unit,  each  of  which  contains  an  automated,  six-position  filter  wheel  to  help  filter  the  excitation  laser 
source  from  the  fluorescence  signal.  The  input  optics  collimate  the  light  from  the  fiber,  pass  it  through  a  filter 
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of  the  user’s  choice,  then  focus  the  light  on  the  spectrometer  entrance  slit  matching  the  f-number  of  the 
spectrometer  system.  All  spectrometers  are  controlled  simultaneously  with  programs  written  in  Labview  which 
cycle  the  laser  source  sequentially  through  each  fiber  while  the  other  15  fibers  detect  the  emitted  light  from  the 
tissue  volume,  resulting  in  a  total  of  240  data  points  for  a  given  acquisition.  Exposure  times  are  determined 
automatically  during  acquisition  to  optimize  the  SNR. 

Imaging  system  characterization 

In  funding  year  two,  a  comprehensive  characterization  of  the  spectrometer-based  fluorescence  imaging  system 
which  included  methodical  assessments  of  measurement  repeatability  and  system  sensitivity  to  fluorescence 

activity  was  completed.  Most  of  the 
information  in  this  section  is  presented  in 
a  recent  publication[l]  (attached  as 
Appendix),  though  the  highlights  are 
provided  here.  Repeatability 

measurements  were  made  for  two 
imaging  modes,  one  records  the  intensity 
of  the  excitation  source  while  the  other 
filters  the  excitation  light  and  measures 
the  fluorescence  emission  intensity.  An 
8.6cm  diameter  homogeneous  phantom 
composed  of  silicone,  titanium  dioxide, 
and  India  ink  was  used  to  measure  the 
repeatability  of  transmission  mode 
measurements  using  the  spectroscopy 
system.  Optical  properties  of  the 
phantom  were  0.004  mm'1  and  1.91  mm'1 
for  the  absorption  and  reduced  scattering 
coefficients  respectively.  The  average 
and  maximum  standard  errors  are  0.28% 
and  0.37%  for  fibers  remaining  in  contact 
with  the  phantom  and  14.8%  andl6.3% 
for  re-positioned  fibers,  indicating  the 
larges  source  of  error  is  fiber  position. 
However,  this  can  be  addressed,  as 
described  below. 

Determining  measurement  repeatability 
in  fluorescence  mode  is  less 
straightforward  given  the  signal  dependence  on  fluorophore  concentration,  absorption  spectrum,  and  quantum 
yield.  For  this  study,  a  70  mm  diameter  liquid  phantom  containing  DPBS,  1%  intralipid,  India  ink.  The 
fluorophore  used  in  this  study  was  indocayanine  green  (ICG),  perhaps  the  most  commonly  used  fluorophore  for 
fluorescence  tomography  in  phantoms.  A  concentration  of  lOnM  ICG  was  used  to  determine  measurement 
repeatability.  Repeatability  was  calculated  in  two  ways,  one  considered  only  the  raw  data  for  a  given  binned 
pixel  array  and  resulted  in  a  mean  standard  error  of  0.7%  and  maximum  standard  error  of  1.4%  for  the  pixel 
array  at  the  fluorescence  peak.  The  second  measure  was  determined  based  on  the  integrated  intensity  from  the 
full  calibration  and  spectral  fitting  routine,  resulting  in  an  average  standard  error  of  0.6%  and  maximum 
standard  error  of  1.8%.  It  is  clear  that  fiber  positioning  variability  would  dominate  the  data  error  for 


Rotating  source 
coupling  stage 


Source  fibers  bifurcated  into 


in  the  MRI  control  room.  Long  fiber  optics  penetrate  ports  in  the  wall,  extend  into 
the  bore  of  the  magnet,  and  couple  to  the  tissue  with  custom  designed  fiber  arrays 
attached  to  a  standard  breast  RF  coil. 
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fluorescence  measurements  in  this  case.  However,  fluorescence  measurement  data  offers  a  unique  opportunity 
to  account  for  fiber  coupling  variability  by  calibrating  the  fluorescence  measurements  to  the  transmission 
measurements  in  the  same  geometry.  This  provides  inherent  stability  to  systematic  error. 


A  70  mm  diameter  liquid  phantom  containing  DPBS,  1%  intralipid  and  India  ink  was  used  to  investigate  the 
overall  sensitivity  of  the  system  to  Indocyanine  green  (ICG)  fluorescence.  Fluorescence  emission  and 
excitation  transmission  measurements  were  recorded  for  each  phantom,  with  a  maximum  allowed  camera 
integration  time  of  120  seconds  applied  to  the  fluorescence  measurements.  The  typical  data  measured  across  the 
emission  spectrum  are  shown  in  Fig.  2  (a)  for  a  strong  fluorescence  signal  and  Fig.  2  (c)  for  a  weak  signal.  The 
spectral  fitting  procedure  was  used  to  recover  the  true  fluorescence  signal  and  the  non-specific  background 
contributions,  as  shown  in  (b)  and  (d),  respectively.  Additionally,  to  quantitatively  compare  the  spectral  fitting 
procedure  to  more  conventional  means  of  fluorescence  filtering,  the  measured,  un-fitted  spectra  were  integrated 
to  simulate  720  nm  long-pass  filtering. 
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Figure  2.  A  spectral  fitting  routine  is  used  to  decouple  fluorescence  spectra  from  the  measured  signal. 
The  spectral  fitting  algorithm  determines  the  relative  contribution  of  the  basis  spectra  [right  column,  (b) 
and  (d)]  which  best  fits  the  measured  spectrum.  Two  examples  are  presented  above,  both  for  a  70  mm 
diameter  homogeneous  phantom  containing  intralipid,  ink  and  ICG  [500  nM  for  graphs  (a)  and  (b)  and 
lOOpM  for  graphs  (c)  and  (d)].  Recovered  values  of  fluorescence  yield  as  a  function  of  ICG 
concentration  in  a  70  mm  liquid  phantom  using  two  methods  to  process  the  recorded  data  are  shown  in 
(e).  One  method  implements  the  spectral-fitting  technique  to  decouple  background  contamination 
across  the  fluorescence  emission  spectrum  while  the  other  simply  integrates  the  measured  spectrum, 
as  would  be  the  case  if  the  experimental  system  relied  on  long  pass  filtering  alone. 


Values  of  fluorescence  yield  recovered  using  the  two  data  pre-processing  techniques  are  plotted  as  a  function  of 
known  ICG  concentration  in  Fig.  2  (e).  The  linear  fit  shown  in  the  figure  was  computed  using  the  spectrally  fit 
data  with  the  y-intercept  forced  to  zero  and  indicates  a  strong  linear  correlation  of  recovered  fluorescence  yield 
and  ICG  concentration  (R2  =  0.99).  At  concentrations  above  1  nM,  the  fluorescence  yield  values  calculated 
using  data  that  was  filtered  only,  with  no  spectral  fitting,  closely  matches  the  spectrally  fit  results.  Fluorescence 
signals  produced  at  these  fluorophore  concentrations  dominate  the  detected  signal,  a  finding  consistent  with  the 
data  shown  in  Fig.  2  (a-d).  However,  fluorescence  yield  values  recovered  with  spectrally  fit  data  maintain  the 
linear  relationship  at  lower  concentrations  than  those  recovered  using  data  without  the  spectral  separation  of 
background  contamination.  Clearly,  the  two  techniques  diverge  at  1  nM.  At  this  concentration,  the  residuals 
between  the  recovered  value  of  fluorescence  yield  and  the  linear  approximation  are  0.4%  and  48%  for  the 
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spectrally  fit  and  filtering-only  approaches,  respectively.  As  the  fluorophore  concentration  drops  below  1  nM, 
both  techniques  lose  the  consistent  linear  response  observed  at  higher  concentrations,  though  the  filtered-only 
shows  even  less  sensitivity  to  changes  in  fluorophore  concentration.  At  500  pM,  the  residuals  increase  to  57% 
and  192%  for  the  spectrally  fit  and  filtered-only  data  processing  responses,  respectively.  Below  this  level, 
accurate  quantification  of  fluorescence  activity  is  unlikely  in  this  phantom  configuration,  however,  the 
spectrally  fit  data  still  shows  a  stronger  response  to  changes  in  fluorophore  concentration  down  to  10  pM.  The 
residuals  calculated  at  the  lowest  concentration  measured  were  almost  5  times  larger  without  spectral  fitting, 
clearly  indicating  a  more  sensitive,  if  not  accurate,  response  of  the  spectral  pre-processing  technique.  It  should 
be  noted  that  since  fluorophore  quantum  yield  is  not  explicitly  known  in  this  solution,  the  calculated  slope  of 
the  linear  regression  does  not  provide  information  on  the  actual  relationship  between  true  and  recovered 
concentration,  however,  the  linearity  itself  is  a  critical  measure  of  system  performance. 

Imaging  clinically  relevant  volumes 

Funding  year  2006-2007  resulted  in  a  series  of  phantom  images  of  fluorescence  yield  in  cases  where  perfect 
drug  uptake  was  assumed  and  demonstrated  dramatic  improvements  in  imaging  performance  when  MRI 
structural  information  was  used  to  guide  the  fluorescence  recovery  procedure  [2],  These  studies  provided  an 
initial  feasibility  validation,  however,  perfect  uptake  experimental  conditions  are  unrealistic  for  clinical 
applications  since  even  targeted  drugs  are  expected  to  produce  background  signals  in  healthy  tissue  regions.  In 
this  report,  images  of  clinically  relevant  breast-sized  tissue  volumes  were  produced  showing  recovered  images 

of  low  contrast  fluorescence  activity. 

Ninety  two  millimeter  diameter  cylindrical  liquid  phantoms 
composed  of  DPBS,  1%  intralipid,  and  India  ink  were  used 
to  investigate  the  imaging  sensitivity  for  cases  of  imperfect 
drug  uptake.  A  22  mm  diameter  cylindrical  plastic  tube  was 
suspended  in  the  larger  phantom  volume  to  simulate  a  tumor 
region  of  elevated  fluorophore  concentration.  The 
background  liquid  contained  300  nM  of  ICG  while  the  tumor 
region  was  composed  of  the  same  intralipid  solution  but  with 
varying  concentrations  of  the  fluorophore,  to  produce  a  range 
of  tumor  to  background  fluorophore  concentrations  between 
6.6:1  to  1.5:1,  or  560%  to  50%  contrast,  respectively. 
Images  were  produced  using  three  reconstruction  techniques. 
One  technique  assumes  only  the  outer  boundary  of  the 
domain  is  known.  The  other  two  techniques  use  the  internal 
structure  of  the  phantom,  as  would  be  determined  by 

simulateneously  acquired  MR  images,  to  guide  the 

fluorescence  recovery.  This  information  is  applied  either  as 
strict  knowledge  of  the  internal  anatomy,  (“hard  priors”),  or 
as  a  regulatization  filter  matrix  which  is  more  flexible  in 
terms  of  guiding  the  solution  (“soft  priors”)  [2], 

Reconstructed  images  using  all  techniques  are  shown  in  Fig.  3.  It  is  clear  that  in  this  relatively  large  domain 

with  imperfect  drug  uptake,  reconstructing  images  without  the  tissue  morphology  information  produces 

qualitatively  inaccurate  and  misleading  images.  The  tumor  region  is  not  qualitatively  or  quantitatively 
recovered  even  at  the  higher  contrast  levels.  The  spatial  priors  techniques,  on  the  other  hand,  essentially  reduce 
the  imaging  problem  from  one  of  both  localization  and  quantification  to  one  of  quantification  only  since  the 


Figure  3.  Reconstructed  images  of  fluorescence 
yield  for  liquid  phantoms  containing  intralipid,  ink 
and  ICG,  using  no  spatial  priors  (a),  soft  priors  (b) 
and  hard  priors  (c).  In  this  case,  the  simulated 
tumor  region  was  close  to  the  edge.  Columns 
correspond  to  different  tumor-to-background  ICG 
contrasts,  1 .5:1 , 3.3:1 ,  and  6.6:1  from  left  to  right. 
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boundaries  of  the  tumor  are  specified  by 
the  MR  image.  In  this  case,  the  tumor 
region  is  clearly  visible  in  the  image, 
even  down  to  the  lowest  contrasts  of 
50%.  The  improvement  provided  by  the 
spatially  guided  reconstruction  is  even 
more  dramatic  for  cases  of  imperfect 
uptake,  presented  here,  than  for  infinite 
contrast  cases.  This  is  an  important 
result  since  in  vivo  concentrations  will 
likely  be  closer  to  those  used  in  this 
study. 

Despite  the  added  difficulty  of  imaging 
imperfect  drug  uptake,  the  phantoms 
used  above  are  still  relatively  simple 
compared  to  the  heterogeneous  structure 
of  breast  tissue.  Optically,  the  organ 
can  be  considered  as  a  composition  of  a 
small  number  of  optically  significant 
tissue  types;  adipose  tissue,  fibro-glandular  tissue,  and  suspect  regions  (benign  and  malignant  tumors,  cysts). 
In  order  to  investigate  the  imaging  capabilities  of  the  fluorescence  tomography  system  in  heterogeneous 
imaging  volumes,  a  more  complex  phantom  was  created  using  gelatin.  The  91mm  diameter  cylindrical 
phantom  was  composed  of  water,  gelatin,  TiC>2  for  scattering,  and  ICG.  The  resulting  phantom  is  pictured  in 
Fig.  4  which  shows  the  three  layers;  the  outer  layer  represents  adipose  tissue,  the  second  layer  fibroglandular 
tissue,  and  the  small  cylindrical  inclusion  is  the  simulated  tumor  region.  ICG  concentrations  provide  a  tumor  to 
adipose  contrast  of  10  to  1  and  a  tumor  to  fibroglandular  contrast  of  around  3.3  to  1.  Gadolinium  was  also 
added  to  the  second  layer  so  that  each  simulated  tissue  region  would  be  discemable  an  MRI  scan  for  use  as 
spatial  prior  information.  MRI  and  optical  data  were  collected  and  the  MR  images  segmented  to  provide  spatial 
guidance  (hard  and  soft  priors  techniques)  of  the  fluorescence  reconstructions.  Resulting  images  are  shown  in 
Fig.  4  for  the  spatially  guided  and  unguided  reconstructions.  As  is  clear  from  the  images,  spatial  information  is 
critical  to  obtain  accurate  images  in  complex  domains. 

Absorption  tomography  for  imaging  exogenous  contrast 

Another  approach  to  image  exogenous  contrast  uses  the  optical  absorption  characteristics  of  a  contrast  agent 
rather  than  the  emitted  light  from  a  fluorescent  contrast  agent.  A  “spectrally  constrained”  approach  which 
recovers  images  of  endogenous  tissue  chromophore  (oxy  and  deoxy  hemoglobin,  water)  concentrations  and 
scattering  parameters  using  the  chromophores’  known  spectral  absorption  profiles  as  prior  information  has  been 
shown  in  previous  studies  to  accurately  image  the  endogenous  chromophore  concentration  in  human  breast  [3] 
[4]  [5,  6].  Under  the  current  grant,  this  approach  has  been  expanded  to  incorporate  exogenous  absorbing  drugs. 
Simulation  studies  have  been  completed  which  demonstrate  simultaneous  recovery  of  endogenous  and 
exogenous  chromophores  and  scattering  parameters  using  multi-wavelength  spectrally-constrained  absorption 
tomography.  The  target  values  of  endogenous  and  exogenous  Lutex  chromophore  concentrations  depicted  in 
Fig.  5  were  used  to  generate  simulated  noisy  frequency  domain  data  at  ten  discrete  wavelengths.  This  data  was 
used  to  complete  a  spectrally  constrained  reconstruction  to  recover  images  of  the  target  domain.  The  results 
presented  in  Fig.  5  indicate  reasonable  qualitative  and  quantitative  accuracy  for  the  recovery  of  both  the 
endogenous  chromophores  as  well  as  LuTex. 


(a)  (b)  (c) 


(d)  {«>  (f) 


Figure  4.  A  photograph  of  the  phantom  in  the  optical  fiber  array  is  shown  in  (a). 
The  second  layer  contained  gadolinium,  which  provided  MR  image  contrast,  (b),  for 
easy  image  segmentation  using  Mimics  software  (c).  Fluorescence  yield  images  of 
the  three  layer  phantom  reconstructed  using  no  spatial  guidance  (d),  soft-priors  (e), 
and  hard-priors  (f).  All  images  are  plotted  using  the  same  color  scale. 
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Figure  5.  Target  and  recovered  values  of  endogenous  and  exogenous  chromophore  in  a  realistic 
breast  volume.  Images  were  reconstructed  using  a  spectrally  constrained  technique  which  incorporates 
the  spectral  absorption  features  as  prior  knowledge. 


Comparing  absorption  and  fluorescence  imaging 


This  study  addressed  the  question  of  which  optical  property  provides  the  best  opportunity  for  imaging 
exogenous  optical  contrast  directly  by  comparing  the  intensity  perturbation  of  boundary  data  caused  by  the 
presence  of  an  exogenous  absorber  to  that  caused  by  an  exogenous  fluorophore.  A  breast-shaped  domain 
derived  from  an  MRI  scan  provides  the  imaging  test  field  for  this  simulation  study.  The  “background” 
domain  contains  no  tumor  region  and  presents  endogenous  chromophore  contrast  only  between  the  fatty  and 
fibro-glandular  regions,  as  shown  in  Fig.  6.  Also  represented  in  the  study  are  two  different  tumor-to-normal 
drug  uptake  situations,  one  an  idealized  case  with  infinite  contrast  [Fig.  6  (a)],  and  the  other  a  more  realistic 
drug  uptake  case  resulting  in  finite  contrast  (ranging  from  2:1  to  3:1)  in  drug  concentration  [Fig.  6  (b)].  In 
both  cases  exogenous  absorption  and  fluorescence  contrast  was  introduced  as  a  simulated  drug  with 
absorbing  properties  of  Lutex.  Fluorescence  quantum  yield  was  varied  as  an  independent  variable  for  the 
fluorescence  emission  measurements. 
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Figure  6.  Two  test  domains  for  investigating  perturbations  caused  by  absorption  and  fluorescence  exogenous 
contrast.  In  (a),  exogenous  contrast  tumor  specificity  is  assumed  to  be  perfect,  while  in  (b),  the  uptake  profile  is  more 
realistic.  The  source-detector  configuration  is  shown  in  (c). 


The  parameter  of  interest  in  this  study  is  the  relative  perturbation  in  boundary  data  intensity  caused  by  a 
drug-enhanced  tumor.  Data  from  perturbations  caused  only  by  exogenous  absorption  are  termed 
“transmission”  measurements,  for  which  735  nm  was  used  as  the  laser  source  wavelength.  Fluorescence 
data  was  generated  for  excitation  and  emission  wavelengths  of  690  and  761  nm  and  the  fluorescence 
emission  perturbation  was  recorded  as  a  function  of  fluorescence  quantum  yield  for  excitation  filtering 
efficiencies  of  3  OD,  5  OD  and  7  OD.  Perturbation  values  were  calculated  for  three  detector  positions 
relative  to  a  single  source,  as  shown  in  Fig.  6  (c),  for  both  transmission  and  fluorescence  measurements. 
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Figure  7  presents  graphs  of  relative  intensity  perturbation  as  a  function  of  fluorescence  quantum  yield  at 
three  detector  positions  for  infinite  exogenous  agent  contrast,  corresponding  to  the  test  domain  in  Fig.  6  (a). 
The  graphs  include  perturbations  of  transmission  intensity  resulting  only  from  an  increase  in  the 
concentration  of  the  exogenous  agent  as  well  as  fluorescence  emission  intensity  degraded  by  different 
amounts  of  excitation  bleed-through.  The  perturbation  caused  by  absorption  is  unaffected  by  quantum  yield 
while  perturbations  in  fluorescence  signals  are  proportional  to  quantum  yield. 


<u 


—  transmission  amplitude 

-  fluorescence  emission  amplitude  3  OD 

-  fluorescence  emission  amplitude  5  OD 

-  fluorescence  emission  amplitude  7  OD 


quantum  yield 


Figure  7.  Perturbations  in  transmission  amplitude  and  fluorescence  amplitude,  given  different  filtering  efficiencies  and  quantum 
yield  values,  caused  by  a  centrally  located  object  with  perfect  drug  uptake.  Under  these  conditions,  fluorescence  signals  appear  to 
be  more  sensitive  to  the  object. 


It  is  clear  from  Fig.  7  that  if  the  tissue  contains  infinite  tumor-to-background  contrast  of  exogenous  agent 
concentration  the  fluorescence  signal  is  more  sensitive  to  the  presence  of  the  drug  within  the  imaging 
domain.  This  is  consistent  for  all  source-detector  pairs  and  filtering  efficiency  of  5  OD  and  higher  given 
quantum  yields  of  0.0001  and  above.  These  results  substantiate  a  conclusion  favoring  fluorescence  imaging, 
however,  many  contrast  agents,  including  all  optical  agents  currently  approved  for  in  vivo  human  use,  will 
not  provide  infinite  specificity  and  imperfect  uptake  must  also  be  considered. 


Figure  8  provides  results  for  the  case  with  drug  uptake  of  2:1  for  tumor-to-fibro-glandular  tissue  layer  and 
3:1  for  tumor-to-fatty  tissue.  For  this  case,  the  absorption  perturbation  is  shown  to  be  more  significant  for  all 
values  of  quantum  yield  and  filtering  efficiencies  tested.  Significantly,  these  values  represent  the  maximum 
expected  perturbation  for  fluorescent  anomalies  and  are  significantly  lower  than  perturbations  caused  by  the 
absorption  profile  of  the  exogenous  drug. 
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Figure  8.  Perturbations  in  transmission  amplitude  and  fluorescence  amplitude,  given  different  filtering  efficiencies  and  quantum 
yield  values,  caused  by  a  centrally  located  object  with  imperfect  drug  uptake.  For  this  lower  drug  contrast  case,  absorption 
measurements  are  more  sensitive  to  the  object  than  fluorescence  emission,  regardless  of  fluorescence  quantum  yield. 


The  quantum  yield  value  at  which  fluorescence  and  transmission  perturbations  intersect,  i.e.  the  quantum 
yield  threshold,  was  determined  for  drug  contrasts  ranging  from  1.1:1  to  10:1,  and  background 
concentrations  ranging  from  10  nM  to  1000  nM.  The  secant  root  finding  method,  initialized  with  the  output 
of  a  bisection  method  algorithm,  was  used  to  determine  the  quantum  yield  value  which  minimizes  the 
difference  between  fluorescence  and  transmission  perturbations  over  the  range  rj  =  10'6  to  t]  =1  for  each 
combination  of  contrast  and  background  concentration.  The  exogenous  contrast  distribution  was  assumed  to 
be  homogeneous  except  for  the  tumor  region  and  only  detector  #4  was  considered  for  this  experiment. 


Threshold  quantum  yield  values  as  a  function  of  background  drug  concentration  for  a  single  contrast  value 

(5:1  in  this  case)  are  presented  in 
Fig.  9.  These  curves  represent  the 
minimum  value  of  quantum  yield 
required  to  ensure  fluorescence 
measurements  are  more  sensitive 
to  the  tumor  region  than 
transmission  measurements  for  a 
range  of  excitation  filtering 
efficiencies.  In  this  manner, 
regimes  may  be  defined  based  on 
contrast  agent  quantum  yield  and 
expected  background 

concentration.  Regions  above  and 
to  the  left  of  a  given  curve 
represent  a  “fluorescence 


(a) 
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Figure  9.  Threshold  quantum  yield  values  plotted  as  a  function  of  background  drug 
concentration  for  a  5:1  tumor-to-background  contrast  (a).  The  quantum  yield 
threshold  curves  delineate  between  experimental  conditions  which  favor 
fluorescence  and  absorption  measurements  (b). 


sensitive”  regime  while  background/quantum  yield  combinations  below  and  to  the  right  of  each  curve  can  be 
considered  “absorption  sensitive”,  assuming  a  5:1  target-to-background  contrast.  This  is  illustrated  in  Fig.  9 
(b)  for  the  7  OD  filter  results.  Conditions  clearly  favor  transmission  measurements  when  background  drug 
concentration  approaches  150  nM,  regardless  of  filtering  efficiency.  On  the  lower  concentration  side  of  this 
hard  limit,  excitation  filtering  efficiency  has  a  large  impact  on  threshold  quantum  yield  values.  As  filtering 
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efficiency  degrades,  higher  and  higher  quantum  yields  are  required  to  produce  meaningful  fluorescence 
perturbations  at  the  tissue  boundary. 
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LuTex  fluorescence  emission  peak  shift  in  turbid  phantoms 

The  influence  of  tissue  optical  properties  on  the  shape  of  near-infrared  (NIR)  fluorescence  emission  spectra 
propagating  through  multiple  centimeters  of  tissue  was  investigated.  The  results  of  this  study  will  be 

appearing  in  the  Journal  of  Applied 
Physics  [7]  and  the  highlights  are 
summarized  here.  Experimentally 

measured  fluorescence  emission  spectra 
measured  in  6  cm  homogeneous  tissue 
phantoms  shows  dramatic  spectral 
distortion  which  results  in  emission  peak 
shifts  of  up  to  60  nm  in  wavelength. 
Measured  spectral  shapes  are  highly 
dependent  on  the  photon  path-length  and 
the  highly  scattered  photon  field  in  the  NIR 
amplifies  the  wavelength-dependent 
absorption  of  the  fluorescence  spectra. 
Simulations  of  the  peak  propagation  using 
diffusion  modeling  describe  the 
experimental  observations  and  confirm  the 
path-length  dependence  of  fluorescence 
emission,  as  shown  in  Fig.  10.  Spectral 
changes  are  largest  for  longer  path-length 
measurements,  and  this  will  be  most 
dominant  in  human  tomography  studies  in 
the  NIR.  Spectrally  resolved  or  multi¬ 
wavelength  band-pass  measurements  are 
required  to  detect  these  changes,  and  may 
be  essential  to  interpret  such  effects,  which 
would  otherwise  be  attributed  to  erroneous 
intensity  measurement.  This  phenomenon 
is  analogous  to  beam  hardening  in  x-ray 
tomography,  which  can  lead  to  image 
artifacts  without  appropriate  compensation. 
The  peak  shift  toward  longer  wavelengths, 
and  therefore  lower  energy  photons, 
observed  for  NIR  luminescent  signals 
propagating  through  tissue  may  readily  be 
described  as  a  “beam  softening” 
phenomenon.  These  spectral  changes 
should  be  considered  for  fluorescence 
tomographic  imaging  through  several 
centimeters  of  tissue.  Use  of  spectrally 
resolved  detection  allows  quantification  of 
this  change,  and  may  be  the  only  reliable 
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Figure  10.  Diffusion  based  modeling  of  the  fluorescence  peak 
through  a  60  mm  diameter  phantom.  The  dilute  emission 
spectrum  of  LuTex  is  shown  as  dotted  lines,  and  the  calculated 
emission  spectrum  from  diffuse  emission  through  the  region  is 
shown  at  discrete  wavelengths  (black  circles).  Experimental 
measurements  are  shown  in  blue  for  intralipid  liquid  phantoms. 
The  circular  domain  shown  represents  a  cross-section  of  the 
cylindrical  phantom  upon  which  the  intensity  (logarithm)  of  the 
excitation  field  calculated  from  the  diffusion  equation  is  plotted  for 
illustration  purposes. 
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way  to  track  intensity  changes  which  would  otherwise  appear  erroneous. 


Contrast-detail  analysis 

A  contrast-detail  study  using  simulated  data  was  completed  to  determine  the  imaging  limits  of  the 
reconstruction  algorithm[8].  Two  circular  domains,  one  51  and  the  other  86  mm  diameter,  were  used  as  test 
fields  and  object  (simulated  tumor)  locations  near  the  edge  and  center  of  each  domain  were  considered.  The 
contrast-detail  results  for  a  threshold  contrast-to-noise  limit  of  three  or  greater  are  plotted  in  Fig.  1 1  for  two 
test  fields  and  object  positions.  Objects  which  are  recovered  with  greater  than  CNR  =  3  have  contrast-detail 

characteristics  which  are  above  and  to  the 
right  of  the  lines  shown  in  Fig.  11  and  are 
thus  considered  detectable,  while  those  below 
and  to  the  left  are  too  small  or  have  too  little 
contrast  to  be  recovered  with  CNR>3  in  the 
image.  Accordingly,  the  limiting  diameter  for 
an  object  near  the  edge  of  the  field  is 
approximately  1.7  mm  for  both  test  field 
diameters.  There  seems  to  be  little  influence 
of  test  field  diameter  on  CNR  limits  for  high- 
contrast  objects  near  the  edge.  Above  an 
object  contrast  of  about  8,  the  continued 
decrease  in  the  minimum  size  for  CNR  =  3,  is 
very  small  with  increasing  contrast,  indicating 
that  this  size  is  a  fundamental  limit  of  the 
imaging  algorithm  for  this  geometry.  This 
high  contrast,  small  detail  portion  of  the  curve 
is  often  referred  to  as  a  “spatial  resolution 
limited”  regime.  It  does  not  appear  that  the 
fundamental  limits  have  been  reached  for 


Figure  11.  Contrast-detail  curve  showing  the  contrast-to-noise 
ratio  (CNR)  of  3  to  approximate  limits  of  detectable  contrast 
and  diameter  for  two  anomaly  positions  in  51  mm  (dashed 
lines)  and  86  mm  (solid  lines)  diameter  test  fields.  In  this 
analysis,  objects  above  and  to  the  right  of  each  line  are  in  the 
region  where  CNR  >  3  and  detection  is  considered  possible. 


objects  near  the  center  of  either  test  field  for  the  contrast  range  studied  and  expected  to  be  encountered 
experimentally.  Furthermore,  CNR  value  of  objects  near  the  center  is  strongly  influenced  by  the  test  field 
size.  For  a  51  mm  test  object,  the  limiting  diameter  is  approximately  4  mm,  for  the  maximum  object  contrast 
(10).  This  increases  to  approximately  8.5  mm  for  the  larger  test  field.  These  results  indicate  the  dramatic 
decrease  in  sensitivity  for  objects  deeper  in  the  test  field. 


Key  research  accomplishments 

•  Constructed  and  calibrated  experimental  system  and  completed  comprehensive  system  performance 
studies  to  determine  signal-to-noise  and  sensitivity  to  contrast  agent  concentration  characteristics. 

•  Incorporated  tissue  structural  information  obtained  from  MR  images  into  the  image  reconstruction 
algorithms  and  demonstrated  a  dramatic  qualitative  and  quantitative  improvement  in  imaging 
performance  using  numerical  and  phantom  data. 

•  Developed  and  tested  algorithms  for  simultaneous  reconstruction  of  endogenous  and  exogenous 
chromophores  based  on  absorption. 

•  Produced  images  of  fluorescence  yield  in  breast-sized  phantoms  with  imperfect  drug  uptake. 

•  Compared  measurement  sensitivity  to  contrasts  arising  from  optical  absorption  and  fluorescence  signals 
in  simulated  breast  domains  and  defined  regions  for  which  each  signal  provides  more  sensitive 
measurements  based  on  contrast  and  overall  concentration. 

•  Observed  fluorescence  emission  peak  shift  propagating  through  turbid  media  and  developed  software  to 
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model  phenomenon  accurately. 

•  Determined  imaging  limits  of  the  numerical  algorithm  using  contrast-detail  analysis. 

Reportable  outcomes 

This  fellowship  led  to  three  first  peer-reviewed  publications,  three  oral  presentations  and  two  poster 

presentations.  These  are  listed  below  along  with  non-first  author  publications  and  conference  proceedings. 

Peer-reviewed  publications 

1.  Dehghani  H,  Davis  SC,  Jiang  SD,  Pogue  BW,  Paulsen  KD,  Patterson  MS,  “Spectrally  resolved 
bioluminescence  optical  tomography,”  Optics  Letters,  3 1(3):365— 367  (2006). 

2.  Pogue  BW,  Davis  SC,  Song  XM,  Brooksby  BA,  Dehghani  H,  Paulsen  KD,  “Image  analysis  methods 
for  diffuse  optical  tomography,”  J.  Biomed  Optics,  99(3):033001 : 1-12  (2006). 

3.  Kepshire  DS,  Davis  SC,  Dehghani  H,  Paulsen  KD,  Pogue  BW,  “Sub-Surface  Diffuse  Optical 
Tomography  can  Localize  Absorber  and  Fluorescent  Objects  but  Recovered  Image  Sensitivity  in 
Non-Linear  with  Depth,”  Applied  Optics,  46(10):  1669-1678  (2007). 

4.  Davis  SC,  Dehghani  H,  Wang  J,  Jiang  S,  Pogue  BW,  Paulsen  KD,  “Image-guided  diffuse  optical 
fluorescence  tomography  implemented  with  Laplacian-type  regularization,”  Optics  Express, 
15(7):4066-4082,  2007. 

5.  Kepshire  D,  Davis  SC,  Dehghani  H,  Paulsen  KD,  Pogue  BW,  “Fluorescence  tomography 
characterization  for  sub-surface  imaging  with  protoporphyrin  IX,”  Optics  Express  16(12):8581-8593 
(2008). 

6.  Davis  SC,  Pogue  BW,  Springett  R,  Leussler  C,  Mazurkewitz  P,  Tuttle  SB,  Gibbs-Strauss  SL,  Jiang 
S,  Dehghani  H, ,  Paulsen  KD,  “Magnetic  resonance-coupled  fluorescence  tomography  scanner  for 
molecular  imaging  of  tissue,”  Review  of  Scientific  Instruments,  79-064302-1-10  (2008). 

7.  Dehghani  H,  Eames  ME,  Yalavarthy  PK,  Davis  SC,  Srinivasan  S,  Carpenter  CM,  Pogue  BW, 

Paulsen  KD,  "Near-infrared  optical  tomography  using  NIRFAST:  Algorithms  for  numerical  model 
and  image  reconstruction  algorithms,”  Communications  in  Numerical  Methods  in  Engineering  , 
(Accepted  2008). 

8.  Wang  J,  Davis  SC,  Srinivasan  S,  Jiang  S,  Pogue  BW,  Paulsen  KD,  “Spectral  tomography  with 
diffuse  near-infrared  light:  inclusion  of  broadband  frequency-domain  spectral  data,”  Journal  of 
Biomedical  Optics  (In  press  2008). 

9.  Davis  SC,  Pogue  BW,  Tuttle  SB,  Dehghani  H,  Paulsen  KD,  “Luminescence  peak  distortion  in 
diffuse  molecular  spectroscopy  of  tissue,”  Journal  of  Applied  Physics,  (Accepted  2008). 

First  author  conference  presentations  and  posters 

1 .  Poster:  “MR-coupled  multi-spectral  diffuse  NIR  imaging  system  for  exogenous  fluorescence  and 
absorption  Tomography,”  Optical  Society  of  America:  Biomedical  Optics  Topical  Meeting,  Ft. 
Lauderdale,  FL  (March,  2006). 

2.  Poster:  “Multi-spectral,  MR-guided  diffuse  optical  tomography  for  imaging  fluorescence  and 
absorption  Exogenous  Contrast,”  Society  for  Molecular  Imaging:  Hilton  Waikoloa  Village,  HI 
(August  2006). 

3.  Presentation:  “MR-guided  reconstruction  of  fluorescence  diffuse  optical  tomography,”  SPIE: 

Photonics  West,  San  Jose,  CA  (January  2007). 

4.  Presentation:  “MR-coupled  molecular  imaging  of  cancer  in  nude  mice  and  breast-sized  domains,” 

SPIE:  Photonics  West,  San  Jose,  CA  (January  2008). 

5.  Presentation:  “MR-coupled  molecular  imaging  of  tissue,”  Optical  Society  of  America:  Biomedical 
Optics  Topical  Meeting,”  St.  Petersburg,  FL  (March,  2008). 
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Conference  Proceedings 

1.  Davis  SC,  Pogue  BW,  Srinivasan  S,  Dehghani  H  and  Paulsen  KD,  "Development  of  spectrally- 
constrained  diffuse  optical  tomography  for  imaging  exogenous  contrast  agents",  Biomedical  Optics 
Topical  Meeting  on  CD-ROM  (The  Optical  Society  of  America,  Washington,  DC),  March  19-22, 
(2006). 

2.  Jiang  S,  Pogue  BW,  Davis  SC,  Paulsen  KD,  "Multispectral  NIR  diffuse  optical  tomography  system 
development",  Biomedical  Optics  Topical  Meeting  on  CD-ROM  (The  Optical  Society  of  America, 
Washington,  DC),  March  19-22,  (2006). 

3.  Kepshire  DS,  Gibbs  SL,  Davis  SC,  Dehghani  H,  Paulsen  KD,  Pogue  BW,  "Diffuse  fluorescence 
tomography  analysis  of  B-scan  mode  geometry",  Biomedical  Optics  Topical  Meeting  on  CD-ROM  (The 
Optical  Society  of  America,  Washington,  DC),  March  19-22,  (2006). 

4.  Carpenter  CM,  Pogue  BW,  Yalavarthy  PK,  Davis  SC,  Jiang  S,  Dehghani  H,  Paulsen  KD,  "Analysis  of 
3-dimensional  reconstruction  in  a  MR-guided  NIR  tomography  system",  Biomedical  Optics  Topical 
Meeting  on  CD-ROM  (The  Optical  Society  of  America,  Washington,  DC),  March  19-22,  (2006). 
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Conclusions 

The  work  completed  as  part  of  this  fellowship  has  produced  a  working  MRI-FMT  system  and  comprehensive 
assessments  of  the  critical  parameters.  Results  using  this  system  demonstrate  a  linear  response  to  fluorophore 
concentration  down  to  1  nM  in  70  mm  cylindrical  phantoms.  Importantly,  images  were  produced  of  complex 
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phantoms  simulating  relatively  low  tumor  to  normal  tissue  contrasts,  i.e.,  poor  drug  specificity,  using  spatially 
guided  reconstruction  techniques  available  only  to  systems  incorporated  into  clinical  imaging  modalities  which 
provide  morphological  information  of  the  breast  tissue.  In  all  cases  studied,  the  multi-modal  imaging  approach 
improved  fluorescence  imaging  performance  dramatically  and  in  many  cases,  was  a  necessary  condition  for 
producing  extracting  any  fluorescence  activity  contrast  in  the  image. 

Simulation  studies  of  realistic  breast  domains  compared  perturbations  caused  by  changes  in  fluorescence  and 
absorption  signals  arising  from  tumor  uptake  of  a  fluorescent  or  absorbing  optical  agent.  Results  indicate  drug 
concentration  is  a  major  factor  determining  which  optical  signal  provides  the  stronger  data  response.  Finally, 
an  observed  fluorescence  peak  shift  in  turbid  media  was  modeled  numerically.  These  spectral  changes  should 
be  considered  for  fluorescence  tomographic  imaging  through  several  centimeters  of  tissue.  Use  of  spectrally 
resolved  detection  allows  quantification  of  this  change,  and  may  be  the  only  reliable  way  to  track  intensity 
changes  which  would  otherwise  appear  erroneous. 
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Abstract:  A  promising  method  to  incorporate  tissue  structural  information 
into  the  reconstruction  of  diffusion-based  fluorescence  imaging  is 
introduced.  The  method  regularizes  the  inversion  problem  with  a 
Laplacian-type  matrix,  which  inherently  smoothes  pre-defined  tissue,  but 
allows  discontinuities  between  adjacent  regions.  The  technique  is  most 
appropriately  used  when  fluorescence  tomography  is  combined  with 
structural  imaging  systems.  Phantom  and  simulation  studies  were  used  to 
illustrate  significant  improvements  in  quantitative  imaging  and  linearity  of 
response  with  the  new  algorithm.  Images  of  an  inclusion  containing  the 
fluorophore  Lutetium  Texaphyrin  (Lutex)  embedded  in  a  cylindrical 
phantom  are  more  accurate  than  in  situations  where  no  structural 
information  is  available,  and  edge  artifacts  which  are  normally  prevalent 
were  almost  entirely  suppressed.  Most  importantly,  spatial  priors  provided 
a  higher  degree  of  sensitivity  and  accuracy  to  fluorophore  concentration, 
though  both  techniques  suffer  from  image  bias  caused  by  excitation  signal 
leakage.  The  use  of  spatial  priors  becomes  essential  for  accurate  recovery  of 
fluorophore  distributions  in  complex  tissue  volumes.  Simulation  studies 
revealed  an  inability  of  the  “no-priors”  imaging  algorithm  to  recover  Lutex 
fluorescence  yield  in  domains  derived  from  T1  weighted  images  of  a  human 
breast.  The  same  domains  were  reconstructed  accurately  to  within  75%  of 
the  true  values  using  prior  knowledge  of  the  internal  tissue  structure.  This 
algorithmic  approach  will  be  implemented  in  an  MR-coupled  fluorescence 
spectroscopic  tomography  system,  using  the  MR  images  for  the  structural 
template  and  the  fluorescence  data  for  region  quantification. 

©  2007  Optical  Society  of  America 
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1.  Introduction 

Imaging  the  spatial  distribution  of  fluorescence  activity  at  depth  in  tissue  is  a  challenging 
problem  that  has  yet  to  impact  clinical  practice.  The  most  popular  approach  is  diffuse  optical 
fluorescence  tomography  (DOFT),  a  model-based  method  which  approximates  photon 
propagation  as  a  diffuse  field  and  computationally  matches  model  parameters  to  measured 
boundary  data.  The  theoretical  framework  is  well  established  [1-5]  and  feasibility  studies 
have  demonstrated  fluorescence  yield  imaging  in  various  phantom  geometries  [6-9] .  Much  of 
the  most  recent  research  has  focused  on  high-throughput,  full  body  imaging  of  small  animals 
[10,  11]  but  no  human  images  have  been  published  to  date.  DOFT  produces  low  resolution 
images  compared  to  standard  clinical  modalities  such  as  x-ray  CT  and  MRI  due  to  the  highly 
scattered  photon  fields  and  relatively  sparse  measurement  sampling  of  the  tissue  volume.  The 
reconstruction  problem  is  ill-posed  and  underdetermined  and  large  heterogeneity  in  tissue 
optical  properties  caused  by  complex  tissue  morphology  further  challenges  the  imaging 
algorithm.  Sensitivity  drops  with  increasing  depth  resulting  in  a  non-linear  responsivity 
across  the  imaging  field,  making  it  more  difficult  to  image  larger  tissue  volumes. 
Experiments  in  larger  volumes  that  more  closely  resemble  a  human  breast  typically  consider 
unrealistically  high  fluorophore  contrasts  and  simple  geometries  [9,  12].  Though  these  studies 
are  important  for  advancing  the  understanding  of  the  modality,  improving  image  resolution 
and  contrast  sensitivity  is  critical  for  identifying  a  clinical  role  for  DOFT. 

DOFT  is  an  extension  of  the  more  widely  studied  diffuse  optical  tomography  (DOT) 
which  suffers  similarly  from  low  resolution  and  depth-dependent  contrast  sensitivity  in  the 
absence  of  spatial  or  spectral  prior  information  to  guide  the  solution.  However,  methods  to 
incorporate  highly-resolved  anatomical  data  obtained  from  standard  clinical  modalities  have 
improved  the  quantification  of  the  optically  derived  images  [13-15].  These  hybrid  approaches 
lead  to  a  conceptually  new  application  of  optical  tomography,  one  in  which  the  highly 
resolved  imaging  system  provides  a  structural  template  upon  which  volumetric  optical 
spectroscopic  images  are  constructed.  This  framework  may  be  applied  to  either  absorption 
and  scatter  spectroscopy,  or  fluorescence  spectroscopy.  Additional  challenges  specific  to  the 
fluorescence  case  include  lower  signal  intensity,  excitation  source  contamination  of  the 
fluorescence  emission  measurements  due  to  filter  leakage,  and  tissue  optical  property  effects 
on  both  the  excitation  and  emission  photon  propagation.  Applying  spatial  guidance 
techniques  to  the  more  complicated  DOFT  problem  may  yield  even  larger  gains  in  imaging 
capability. 

This  paper  introduces  a  method  to  dramatically  improve  fluorescence  imaging  at  depth  in 
tissue  by  coupling  MRI  and  fluorescence  tomography.  Tissue  structural  information 
determined  from  standard  T1  and  T2  MR  images  is  encoded  as  a  spatial  filter  in  the  DOFT 
reconstruction  algorithm  and  used  to  guide  the  recovery  of  fluorescence  activity.  In  this 
implementation,  reconstruction  parameters  are  loosely  grouped  into  regions  based  on  tissue- 
type  determined  from  the  MR  images  but  are  permitted  to  update  independently,  giving  rise  to 
the  term  “soft”  spatial  priors  [16].  The  algorithm  is  tested  with  phantom  data  of  Lutetium 
Texaphyrin,  a  near-infrared  photosensitizer  shown  to  accumulate  preferential  in  a  variety  of 
malignancies  [17-20],  recorded  from  a  newly  developed  MR-coupled  spectroscopy-based 
fluorescence  tomography  system.  A  simulation  study  based  on  realistic  geometries  generated 
from  MR  images  of  a  normal  human  breast  serves  as  an  initial  illustration  of  expected 
improvements  provided  by  the  spatial  priors  approach  in  complex  domains. 
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2.  Theory 

Assuming  highly  diffuse  media,  the  transport  of  light  at  a  given  modulation  frequency  to  in 
the  presence  of  fluorescence  generated  by  an  external  source  at  the  excitation  wavelength  (^x) 
of  the  fluorescing  agent,  is  modeled  by  two  diffusion  equations,  where  the  solution  of  the  first 
equation  provides  the  driving  source  term  of  the  second  [5]: 

-V-rt(r)VOi(r,ft))  +  (//fl  (r)  +  -^-)Ot(r,«)  =  q0(r,(o)  (1) 

c(r) 

-  V  ■  Km (r)VOm  (r,  <y)  +  (ft  (r)  +  (r,  a)  =  O  c(r,  co)tj<u  (r)  (2> 

m  c(r )  f  1  +  [ C0T(r)] 

where  subscripts  x  and  m  represent  the  excitation  and  emission  fluence  at  wavelengths  Xx  and 
Xm,  respectively.  The  intrinsic  optical  parameters  jda  and  jil's  are  the  absorption  and 

reduced  scattering  coefficients  respectively,  q0(r,CD)  is  an  isotropic  source  and 
is  the  photon  fluence  rate  at  position  r.  The  diffusion  coefficient  is  given  by 

1 

Kx  m  — - ; -  and  c(r)  is  the  speed  of  light  in  the  medium  at  any  point,  defined 

by  Co/n(r ),  where  n(r)  is  the  index  of  refraction  at  the  same  location  and  cQ  is  the  speed  of  light 
in  vacuum.  The  fluorescence  parameters  are  the  lifetime  T(r)  and  the  fluorescence  yield 

TjJLlaf  (r) ,  the  latter  a  product  of  the  fluorophore’s  quantum  efficiency  rj  and  its  absorption 
coefficient  jla  (r) . 

The  most  appropriate  description  of  the  air-tissue  boundary  is  derived  with  an  index- 
mismatched  type  III  condition  (also  known  as  Robin  or  mixed),  in  which  some  fraction  of  the 
fluence  at  the  external  boundary  of  the  tissue  exits  and  does  not  return.  The  flux  leaving  the 
external  boundary  is  equal  to  the  fluence  rate  at  the  boundary  weighted  by  a  factor  that 
accounts  for  the  internal  reflection  of  light  back  into  the  tissue.  This  relationship  is  described 
in  the  following  equation: 

<b,JZ,c»)+2An  =  0  (3) 

where  ^  is  a  point  on  the  external  boundary,  and  A  depends  upon  the  relative  refractive  index 
(RI)  mismatch  between  tissue  Q  and  air.  A  can  be  derived  from  Fresnel’s  law: 

2/(l-i?n)-l  +  Icos  6r  1 3 

A  =  ^ - oj - L - (4) 

1  -  |cos  6C  | 

where  6 c  —  arc  sin  (nAIR  /nx )  ,  the  angle  at  which  total  internal  reflection  occurs  for  photons 

„  (nx  /  nAIR  - 1)2 

moving  from  region  £2  with  RI  nx  to  air  with  RI  tiair,  and  /vn  =  - - —  .  At  the 

(», /«,«+!) 

external  boundaries,  nAIR=  1. 

2. 1  Finite  element  implementation: 

When  the  refractive  index  is  homogenous  (assumed  to  be  1.33  [21]),  the  finite  element 
discretization  of  a  volume  Q  can  be  obtained  by  subdividing  the  domain  into  D  elements 
joined  at  V  vertex  nodes.  In  the  finite  element  method  (FEM),  m  (r)  is  approximated  by 
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the  piecewise  continuous  polynomial  function  Q>hxm(r,w)  =  ,  where 

Q/7  is  a  finite  dimensional  subspace  spanned  by  basis  functions  {ui  (r);  i  =  1 . .  .V }  chosen 

/j 

to  have  limited  support.  The  problem  of  solving  for  <Pxm  becomes  one  of  sparse  matrix 

inversion,  and  in  this  work,  a  bi-conjugate  gradient  stabilized  solver  is  used.  As  developed 
previously  [22,  23],  the  coupled  diffusion  Eqs.  (1)  and  (2)  in  the  FEM  framework  can  be 
expressed  as  a  system  of  linear  algebraic  equations: 

f  -•  ~ 


*,(*•)+ c, <a„  +j^)> +  i®, = e„ 

f  •  i  \ 

®w=-a 

zA  J 


(5) 

(6) 


ICO 


where  the  matrices  Kxm(ic),  Cxm  (//fl  H - ) ,  and  Fxm  have  entries  given  by: 

c(r) 

Kx,mij  -  J  Kx,m(r)Vui(r).Vuj(r)dnr 

Q 

Q,m  =  j (Jla xm  {r)  +  -^-)ui{r)uj(r)dnr 
,  j  «  C(r) 

=  §ui{f)uj(r)dn-lr 

and  the  source  vector  Q0  has  terms 

Q0  =  jM;(r)^r0(r)J"r 


(7) 

(8) 

(9) 

(10) 


The  source  term  is  defined  as  a  Gaussian  distribution,  matching  the  intensity  profile  at  the  tip 
of  the  optical  fiber.  Because  the  source  is  assumed  spherically  isotropic,  modeling  is  more 
accurate  when  it  is  centered  one  scattering  distance  within  the  outer  boundary.  The  source 
vector  Qm  for  fluorescence  re-emission  is  expressed  as 


Qm ,  =ju,(r) 


O  (r  (0)7111  (r)- — ^22 


\dnr 


(11) 


and  is  distributed  throughout  the  domain. 
2.2  The  inverse  model 


2.2.1  No  spatial  priors 

In  the  inverse  problem,  the  goal  is  the  recovery  of  optical  properties  at  each  FEM  node  using 
surface  measurements  of  light  fluence  at  both  the  excitation  and  emission  wavelength 
sequentially  (assuming  the  use  of  two  externally  applied  sources,  one  at  each  wavelength), 
followed  by  fluorescence  yield  reconstruction  at  the  emission  wavelength.  The  computational 

approach  is  to  minimize  the  difference  between  measured  fluence,  ,  at  the  tissue 

surface  and  calculated  data,  <s>cxm  ,  from  the  model  Eqs.  (1)  and  (2)  by  adjusting  the  spatial 

distribution  of  the  unknown  parameters  through  minimization  of  the  ‘objective’  function.  The 
objective  function  for  recovering  the  optical  properties  at  the  excitation 

wavelength,  jlx  —  ( jUa  ,  ) ,  is  given  as 
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NM  .  .  NN  /  \ 

i=l  7=1 


(12) 


where  NM  is  the  total  number  of  measurements  given  by  the  imaging  system,  NN  is  the 
number  of  parameters  representing  the  optical  property  distribution  which  corresponds  to  the 
number  of  nodes  in  the  reconstruction  mesh,  and  I  is  an  NN  x  NN  identity  matrix.  In  general, 

2  ?>X2 

X  will  not  equal  zero,  but  the  values  of  jlx  for  which  -  is  close  to  zero  can  be 

determined,  based  on  an  initial  estimate  of  [lx  .  Following  the  Taylor  series  method  for 


deriving  Newton’s  method, 


dx 2_ 


is  evaluated  at  jUx  based  on  an  expansion  around  some 


nearby  point  JUXq  ,  where  the  second  and  higher  order  terms  are  ignored,  leading  to  the 
iterative  update  equation: 

AjUx  =[jT J  +  hi]' JT (<&Mxeas  -O^)  (13) 

where  J  is  the  Jacobian  matrix,  here  calculated  using  the  Adjoint  method  [24].  Equation  (13) 
is  known  as  the  Moore-Penrose  generalized  inverse  and  is  found  to  be  suitable  for  problems 
where  the  number  of  unknowns  to  be  recovered  is  much  larger  than  the  amount  of  information 
(measurements)  available.  In  standard  practice,  I  is  an  identity  matrix,  and  in  this  work  X  is 
some  fixed  fraction  multiplied  by  the  maximum  value  on  the  diagonal  of  the  Hessian  matrix 

JT  J  ,  and  is  therefore  updated  at  each  iteration.  To  recover  the  optical  properties  at  the 
emission  wavelength,  JLlm  =  [jUa  ,  Km  ) ,  the  externally  applied  illuminating  source  is  changed 

to  one  at  the  emission  wavelength  and  the  formulation  presented  in  Eq.  (13)  is  used. 

The  recovered  optical  properties  are  used  in  the  fluorescence  yield  reconstruction  which 
also  adheres  to  the  minimization  formulation  presented  in  Eq.  (12).  The  unknown  parameter 

in  Eq.  (13)  becomes  7) '/LLa  (r) ,  and  the  Jacobian  can  be  calculated  by  similar  Adjoint 

properties  described  above  and  in  Ref.  [5]. 


2.2.2  Spatial  priors 


Spatial  prior  information  is  incorporated  by  assuming  a  ‘generalized  Tikhonov’  penalty  term 
which  is  similar  in  structure  to  Eq.  (12)  except  that  the  identity  matrix  is  replaced  with  a 
Laplacian-type  matrix,  presented  here  for  the  excitation  field: 

NM  NN  /  x 

z2 =£(*"“-*()  04) 

i=  1  7=1 

where  NN  is  the  number  of  unknowns  in  the  model  [16,  25,  26].  The  constant,  / 3  ,  balances 
the  effect  of  the  parameters  with  the  model-data  mismatch  in  the  same  manner  as  X  in  Eq. 
(13).  The  dimensionless  ‘filter’  matrix,  L,  is  generated  using  MRI-derived  priors  and  its 
construction  is  flexible.  In  this  application,  each  node  in  the  FEM  mesh  is  labeled  according 
to  the  region,  or  tissue  type,  with  which  it  is  associated  (in  the  MR  image).  The  L-matrix 
represents  a  Laplacian-type  structure,  the  diagonal  of  which  is  Lu=  1  where  i  is  the  nodal 
index.  When  nodes  i  and  j  are  in  the  same  region  containing  n  nodes,  L?J=-  1/n,  otherwise 
Lij= 0.  This  effectively  relaxes  the  smoothness  constraints  at  the  interface  between  different 
tissues,  in  directions  normal  to  their  common  boundary.  The  effect  on  image  quality  is  similar 
to  that  achieved  through  total  variation  minimization  schemes  [27]  but  easily  encodes  internal 
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boundary  information  from  MR  images.  Following  the  same  procedure  as  in  the  no-priors 
case,  the  parameter  update  is  given  by 

A nx  =  [jTJ  +  /3Lt l]'1  Jt (of - ®cx )  (15) 

where  much  like  X  in  Eq.  (12),  [3  is  a  fixed  fraction  multiplied  by  the  maximum  value  on  the 

diagonal  of  JT  J  .  This  update  formulation  is  also  used  to  recover  the  optical  properties  at  the 
emission  wavelength  as  well  as  the  fluorescence  yield,  as  described  in  Section  2.2.1. 

2.2.3  Reconstruction  basis 

A  number  of  different  strategies  for  defining  the  reconstruction  basis  are  possible,  including 
second  mesh  and  pixel  basis  [28,  29].  The  choice  of  reconstruction  basis  allows  for 
computational  efficiency  which  serves  to  reduce  the  number  of  unknowns  in  the  algorithm. 
The  problem  at  hand  is  twofold:  The  forward  problem  requires  that  the  volume  of  interest  is 
subdivided  into  adequate  number  of  sub-domains  which  allow  for  an  accurate  description  of 
the  calculated  fields,  whereas  a  reduction  in  the  number  of  unknowns  improves  the  ill- 
posedness  of  the  problem  in  the  reconstruction  algorithm.  This  is  addressed  by  defining  a 
separate  reconstruction  basis  (different  from  the  meshes  used  in  the  FEM  implementation), 
upon  which  the  unknowns  are  updated  in  Eqs.  (13)  and  (15)  [23].  In  cases  where  no  prior 
structural  information  is  available,  a  pixel  basis  which  defines  a  set  of  regularly  spaced  pixels 
for  the  update  of  the  quantities  of  interest  is  used.  However,  when  spatial  priors  are  involved, 
it  is  crucial  to  ensure  that  enough  pixels  are  used  to  adequately  define  small  structural  regions. 
In  this  case,  a  set  of  regular  pixel  bases  are  introduced  in  each  region  of  interest,  which 
account  for  the  region’s  individual  shape  and  size.  A  semi-adaptive  method  which  allows  the 
number  of  pixels  in  each  region  to  be  selected  was  used  in  the  studies  presented  here.  The 
same  region-based  pixel  basis  is  used  throughout  the  iterative  reconstruction  process. 

3.  Methods 

3. 1  Simulation  studies 

Test  domains  for  simulation  studies  were  derived  from  a  T1  weighted  MR  image  of  a  human 
breast  which  measured  approximately  10.5  cm  in  diameter.  A  2-D  slice  of  the  breast  volume 
was  discretized  into  a  mesh  of  approximately  2000  nodes  and  MR  image  intensity  thresholds 
were  used  to  assign  adipose  and  fibro-glandular  tissue  volumes  into  mesh  regions.  Figure  1 
shows  the  original  MR  breast  image  and  its  associated  discretized,  region-labeled  mesh.  The 
MR  image  originated  from  a  clinical  exam  with  our  MR  coupled  frequency  domain  NIR 
system  and  therefore  reveals  an  irregularly  shaped  breast  boundary  caused  by  the  fiber  optic 
array  slightly  compressing  the  breast  at  the  contact  positions.  A  cancerous  tumor  region  was 
added  as  a  target  anomaly  or  region  of  interest  for  these  studies  and  is  depicted  in  the  figure. 
A  second  test  case  with  a  tumor  region  located  near  the  center  of  the  domain  was  used  to 
explore  the  nonlinear  sensitivity  of  diffuse  optical  tomographic  techniques.  The  source- 
detector  positions  were  determined  directly  from  the  MR  image  and  represent  16  optical  fibers 
circumscribing  the  breast  domain  in  a  single  plane.  Light  is  detected  at  all  non-source  fiber 
positions  providing  a  total  of  240  measurements  of  intensity  and  phase  per  wavelength.  This 
simulated  configuration  matches  our  experimental  fluorescence  tomography  system. 
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Fig.  1.  An  axial  T1  weighted  MR  slice  of  a  human  breast  is  shown  in  (a),  from  which  the  test 
domain  for  simulation  studies  was  derived.  Darker  regions  indicate  fibro- glandular  tissue 
imbedded  in  adipose  tissue  indicated  by  lighter  values.  The  image  was  acquired  during  a 
clinical  exam  with  one  of  our  MR- coupled  NIR  tomography  systems  and  shows  the 
indentations  caused  by  the  fiber  optic  probes.  The  test  domain  (b)  was  a  discretization  of  the 
MR  image  thresholded  into  regions.  The  yellow  anomaly  was  added  to  simulate  a  targeted 
cancerous  tumor. 


In  this  set  of  studies,  regions  were  assigned  values  in  terms  of  biologically  relevant  optical 
absorbing  (HbO,  dHb,  Water,  Lutetium  Texaphyrin)  as  well  as  scattering  parameters  (Mie 
scattering  amplitude  and  power).  These  values,  provided  in  Table  1,  represent  typical  known 
in  vivo  levels  of  endogenous  chromophores  and  are  consistent  with  previous  clinical  work 
[30,  31].  Though  it  is  unknown  exactly  how  Lutetium  Texaphyrin  (LuTex)  will  distribute  in  a 
human  breast,  the  concentrations  used  are  similar  to  those  reported  in  ex  vivo  studies  [17]  and 
at  least  represent  a  reasonably  complex  distribution  for  demonstrating  recovery  of 
fluorescence  yield.  Tissue  chromophore  concentrations  were  used  to  calculate  total  optical 
absorption  and  scattering  values  at  the  excitation  and  emission  wavelengths  based  on 
experimentally  determined  values  of  molar  extinction  coefficients.  Simulated  noisy  data  was 
generated  based  on  these  optical  properties  in  the  following  manner:  1)  frequency  domain 
data  for  a  light  source  at  the  excitation  wavelength,  2)  frequency  domain  data  for  a  light 
source  at  the  emission  wavelength,  and,  3)  Continuous  Wave  (CW)  fluorescence  emission 
data  at  the  emission  wavelength.  This  represents  a  total  of  3  data  sets  for  a  given  imaging 
session:  two  frequency  domain  data  sets  for  determining  background  optical  properties  and 
one  CW  fluorescence  emission  data  set.  CW  fluorescence  emission  data  was  used  to  match 
the  capabilities  of  our  experimental  system  and  results  in  a  simplification  of  Eq.  (2)  which  can 
be  handled  by  setting  co  =  0. 


Table  1 .  Chromophore  concentrations  and  scattering  parameter  values  assigned  to  the  mesh  regions  in  the  simulation 
studies 


Oxy¬ 

hemoglobin 

(nM) 

Deoxy¬ 

hemoglobin 

(HM) 

Water 

(%) 

Lutetium 

Texaphyrin 

(pM) 

Scattering 

Amplitude 

Scattering 

Power 

Adipose 

10 

10 

50 

0.3 

1.0 

1.0 

Fibro- 

15 

15 

60 

0.5 

1.1 

1.1 

glandular 

Tumor 

20 

18 

90 

1.0 

1.2 

1.15 

Given  the  modest  Stoke’ s  shift  of  LuTex,  depicted  in  Fig.  2,  the  choice  of  excitation 
wavelength  has  practical  implications  for  experimental  work.  The  difficulty  in  filtering  the 
excitation  signal  precludes  the  use  of  Lutex’s  NIR  absorption  peak  (about  735  nm)  to  excite 
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the  fluorophore.  In  this  study,  a  690  nm  excitation  wavelength  was  used  for  both  simulated 
and  experimental  data  acquisition.  In  addition  to  normally  distributed  random  noise  added  to 
the  frequency  domain  data  (5%  amplitude  and  1  degree  phase)  and  CW  fluorescence  emission 
(10%  intensity),  excitation  signal  leakage  through  the  filter  was  added  to  the  CW  fluorescence 
emission  intensity  to  simulate  a  typical  7  OD  rejection  of  the  excitation  intensity.  This 
number  comes  directly  from  experimentally  measured  rejection  estimates  for  the  filters  in  the 
tomography  system  used  in  the  experiments  performed  here. 


Fig.  2.  The  normalized  absorbance  and  fluorescence  emission  spectra  of  Lutetium  Texaphyrin 
are  shown. 


The  general  image  reconstruction  protocol  was  as  follows: 

1)  Reconstruct  for  optical  properties  at  the  excitation  wavelength,  jaax  and  jlisx’,  with 
frequency  domain  data, 

2)  Reconstruct  for  optical  properties  at  the  emission  wavelength,  jaam  and  jasm’,  with 
frequency  domain  data  collected  using  a  laser  source  at  the  emission  wavelength, 

3)  Use  the  reconstructed  optical  properties  and  fluorescence  intensity  data  to  recover 
fluorescence  yield. 

The  same  reconstruction  algorithm  was  used  to  determine  background  optical  properties  in 
steps  (1)  and  (2)  and  is  based  on  previously  reported  work  [32,  33].  Initial  estimates  for  all 
parameters  were  generated  using  homogenous  fitting  algorithms  which  enforce  a  single  value 
for  all  nodes.  The  Jacobian  matrix  was  calculated  on  a  fine  mesh  of  approximately  2000 
nodes  and  interpolated  onto  a  course  reconstruction  pixel  basis  for  inversion.  A  30  by  30 
pixel  reconstruction  basis  was  used  for  the  no-priors  case.  The  spatial  priors  reconstruction 
used  a  newly  developed  semi-adaptive  pixel  basis  that  redistributes  the  pixels  based  on  the 
region  information,  as  described  in  section  2.2.3.  This  method  ensures  that  each  region 
contains  an  adequate  number  of  nodes  to  approximate  the  internal  structure  of  the  domain. 
Convergence  was  defined  as  less  than  a  2%  change  in  projection  error  between  successive 
iterations  for  the  frequency  domain  optical  properties  algorithm  and  less  than  a  1%  change  in 
projection  error  between  iterations  for  the  fluorescence  yield  reconstruction.  Similar 
algorithmic  parameters  were  used  for  the  phantom  studies  and  are  described  in  further  detail 
below. 

3.2  Phantom  studies 

A  spectrometer  based  tomographic  imaging  system  which  couples  directly  into  a  Philips  3T 
MRI  magnet  was  used  to  acquire  fluorescence  emission  spectra.  The  system,  depicted  in  Fig. 
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3,  is  composed  of  16  spectrometers  with  low  noise  CCD  cameras  cooled  to  -70  C.  Sixteen  13 
meter  long  silica  fiber  bundles  circumscribe  the  imaging  domain  in  contact  mode  and  couple 
into  the  spectrometers  via  custom  designed  input  optics  with  a  six  position  automated  filter 
wheel.  Fluorescence  emission  signals  are  processed  with  long  pass  interference  filters  before 
entering  the  spectrometer.  Each  fiber  bundle  contains  eight  400  jam  fibers,  seven  of  which 
directly  couple  the  tissue  surface  to  the  spectrometer  input  while  the  eighth  couples  the  tissue 
surface  to  the  light  source.  This  detector  configuration  minimizes  coupling  losses  and 
provides  parallel  detection  of  full  spectra  for  each  source  position.  Inter-fiber  tissue  coupling 
calibration  factors  were  determined  prior  to  imaging  using  a  cylindrical  homogenous  phantom 
with  a  centrally  located  source/detector.  The  custom  manufactured  input  optics  were  designed 
to  optimize  filtering  and  maximize  the  light  detection  efficiency  of  the  spectrometers.  Camera 
exposure  times  can  be  adjusted  to  the  detected  light  level  to  maximize  the  signal  to  noise  ratio. 
With  15  detectors  per  16  source  positions,  a  total  of  240  measurements  are  recorded  for  a 
given  acquisition. 


Fig.  3.  The  experimental  spectrometer-based  system  depicted  at  left  couples  directly  into  the 
MR  via  13  meter  fiber  optic  bundles.  Sixteen  spectrometers  are  computer  controlled  for  rapid 
image  acquisition  (left  photograph).  An  animal  interface  (right  photograph)  is  composed  of  a 
rodent  MR  coil  custom  built  by  Philips  Research  Hamburg  to  accommodate  the  optical  fiber 
array  for  simultaneous  MR  and  NIR  fluorescence  imaging. 


Lutetium  Texaphyrin  provided  by  Pharamcyclics  was  diluted  in  water  and  used  as  the 
imaging  fluorophore.  The  test  domain  was  a  solid  5.5  cm  diameter  hardened  epoxy  phantom 
with  scatterer  and  absorbers  created  by  titanium  dioxide  powder  and  India  ink,  respectively 
[34,  35].  The  phantom  had  a  14  mm  hole  located  approximately  12  mm  from  the  phantom 
center.  The  background  optical  properties  were  juax  =  0.005  and  jlisx’  =  1.0  mm'1,  measured 
with  a  frequency  domain  system  near  the  excitation  wavelength.  Unlike  the  simulation 
experiments,  the  optical  properties  were  assumed  constant  throughout  the  domain  in  this 
experiment.  These  values  were  also  used  as  the  optical  properties  at  the  emission  wavelength. 
The  hole  was  filled  with  a  solution  of  1%  Intralipid  to  match  the  scattering  value  of  the 
background,  and  varying  concentrations  of  Lutex  (0.3125  jliM  to  5  juM)  were  added.  This 
represents  a  simple  test  case  for  investigating  the  imaging  response  to  varying  concentrations 
of  fluorophore.  The  excitation  source  was  a  690  nm  laser  diode  which  matches  the 
wavelength  used  in  the  simulation  studies.  Total  acquisition  time  for  the  fluorescence 
emission  was  less  than  4  minutes  (total  of  240  data  points). 

Even  after  processing  the  collected  light  with  a  720  nm  long  pass  interference  filter 
(Omega  Optics)  which  provides  7  OD  rejection  of  the  excitation  light  as  well  as  the  filtering 
offered  by  the  spectrograph  grating,  emission  spectra  recorded  by  the  detector  are  composed 
of  a  sum  of  the  pure  fluorescence  signal  and  excitation  bleed-through.  In  order  to  decouple 
these  signals,  previously  recorded  “basis  spectra”  of  the  bleed-through  signal  and  the  pure 
fluorescence  signal  are  fit  to  the  data.  The  process  is  illustrated  in  Fig.  4  for  one  measurement 
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at  a  single  detector.  A  linear  least  squares  algorithm  operating  on  the  basis  spectra  quantifies 
the  amount  of  excitation  bleed-through  and  true  fluorescence  signal  that  exists  in  each  spectral 
recording,  further  reducing  the  effect  of  excitation  bleed-through.  This  is  accomplished  by 
minimizing  the  summation 

5  =  Xb’, -(«m)  +  ^Ga,))]2  (16) 

i= 1 

with  respect  to  a  and  b,  where  yt  is  the  measured  intensity  at  a  given  wavelength  pixel,  F  and 
G  are  the  excitation  and  fluorescence  basis  spectra,  a  and  b  are  the  coefficients  recovered  in 
the  minimization  procedure,  and  N  is  the  number  of  wavelength  pixels  per  spectrum.  Once 
fit,  the  fluorescence  signal  is  integrated  and  becomes  the  fluorescence  emission  intensity  data 
for  the  reconstruction  algorithm. 


(a) 


(b) 


CCD  Pixel  number 


CCD  Pixel  number 


Fig.  4.  An  example  of  a  pair  of  basis  spectra  for  the  excitation  and  fluorescence  light  (in 
counts/s  as  a  function  of  CCD  pixel  number)  is  shown  in  (b).  These  spectra  are  recorded  for 
each  detector  prior  to  imaging.  In  practice,  the  basis  spectra  are  used  to  perform  a  least  squares 
fit  (a)  to  the  spectrum  measured  for  each  source-detector  pair  to  determine  the  relative 
contribution  of  the  fluorescence  and  excitation  light  to  the  measured  response. 


To  investigate  the  improved  imaging  capability  of  reconstructing  data  with  spatial  priors, 
each  data  set  was  reconstructed  with  and  without  spatial  “soft”  priors.  Since  the  domain  was 
easily  characterized  geometrically,  spatial  prior  information  was  determined  by  direct  manual 
measurement  of  the  phantom.  This  information  was  encoded  in  the  fine  mesh  used  for 
reconstruction.  Since  background  optical  properties  were  previously  determined  with  our 
frequency  domain  system,  they  were  assumed  known  for  image  reconstruction.  Convergence 
was  defined  as  less  than  a  1  %  change  in  projection  error  between  iterations.  All  images  were 
reconstructed  using  a  2GHz  Centrino  Duo  laptop  with  2GB  RAM  running  Windows  XP. 

4.  Results 

4.1  Simulation  results 

Figures  5  and  6  display  the  target  values  of  the  two  test  cases  along  with  images  reconstructed 
using  no  priors  and  spatial  soft  prior  information.  Qualitatively,  image  recovery  using  spatial 
priors  produces  significantly  more  accurate  images.  Spatial  priors  preserve  the  general 
internal  structure  of  the  Lutex  distribution,  detail  that  is  lost  almost  entirely  in  the  no-priors 
case.  Images  of  fluorescence  yield  show  the  most  dramatic  difference  between  the  spatial 
priors  and  no  priors  reconstructions.  Without  spatial  priors,  the  algorithm  appears  to  have  no 
ability  to  recover  “cancer”  regions  of  elevated  fluorescence  yield  for  these  complicated  cases. 
However,  incorporating  spatial  priors  results  in  images  that  qualitatively  appear  accurate  and 
quantitatively  are  reasonably  close  to  the  true  values.  Figure  7  provides  1-D  cross-sections 
near  the  y-axis  of  each  domain.  These  plots  confirm  an  inability  of  the  no-priors  imaging 
algorithm  to  recover  the  simulated  fluorescent  tumor  in  either  test  field.  Alternatively,  the 
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spatial  prior-based  imaging  algorithm  not  only  picks  out  the  objects  of  interest,  but  provides 
fairly  accurate  reproductions  of  the  complicated  structure  of  simulated  fibro-glandular  layers. 
Mean  values  for  the  simulated  cancer  region  near  the  domain  center  are  79%  and  51%  of  the 
true  values  for  the  spatial  priors  and  no-priors  reconstructions,  respectively.  These  numbers 
change  to  75%  and  45%  for  the  case  with  the  cancer  region  closer  to  the  edge  of  the  domain. 
They  represent  a  significant  improvement  in  imaging  performance;  however,  they  alone  do 
not  illustrate  the  full  impact  of  incorporating  spatial  priors.  The  cross-sectional  plots  indicate 
that  the  no-priors  images  contain  virtually  no  spatial  discrimination  of  the  cancer  regions. 
Regional  contrasts  are  depicted  in  Table  2  and  further  illustrate  a  dramatic  overall 
improvement  in  cancer  region  quantification  with  spatial  priors. 


Target  distribution  to  No  spatial  priors  Spatial  “soft”  priors 

generate  data  reconstructed  image  reconstructed  image 


Fig.  5.  Target  and  recovered  values  of  pajX,  p7^,  pajm,  p7Sjm,  and  fluorescence  yield,  r|paf  ,  for 
reconstruction  implementations  using  no  prior  information  and  with  spatial  prior  information. 
In  this  case,  the  simulated  cancer  region  is  near  the  edge,  which  is  known  to  be  easier  to 
recover  without  spatial  priors.  The  image  scales  are  at  right. 
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Target  distribution  to  No  spatial  priors  Spatial  “soft”  priors  j 

generate  data  reconstructed  image  reconstructed  image 


Fig.  6.  Target  and  recovered  values  of  pajX,  p7^,  pajm,  p7Sjm,  and  fluorescence  yield,  r|paf  ,  for 
reconstruction  implementations  using  no  priors  and  spatial  priors.  In  this  case,  the  simulated 
tumor  region  is  near  the  center  of  the  imaging  domain,  which  is  known  to  be  more  difficult  to 
recover  accurately.  Image  scales  are  at  right. 
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Table  2.  Target  and  recovered  fluorescence  yield  regional  contrasts  for  the  images  in  Figs.  4  and  5. 


Shallow  tumor 

Target  Contrast 

No  spatial  priors 

Spatial  “soft”  priors 

Tumor  :  fibro- 

2.0:  1 

1.1  :  1 

1.6  :  1 

glandular 

Tumor  :  adipose 

3.3  :  1 

1.3  :  1 

2.3  :  1 

Deep  tumor 

Target  Contrast 

No  spatial  priors 

Spatial  “soft”  priors 

Tumor  :  fibro- 

2.0:  1 

1.2  :  1 

1.8  :  1 

glandular 

Tumor  :  adipose 

3.3  :  1 

1.4  :  1 

2.4:  1 

Reported  ratios  were  calculated  using  the  mean  values  in  each  region. 

Fig.  7.  Cross  sectional  plots  of  fluorescence  yield  are  shown  for  the  simulated  imaging 
domains  in  (a)  the  case  with  an  object  near  the  edge  and  (b)  the  case  with  an  object  near  the 
center.  In  both  cases,  the  cross  section  is  in  the  y-direction  just  off  center  from  x  =  0.  The 
solid  line  represents  the  target  value,  the  small  dotted  line  the  recovered  value  using  a  no-priors 
based  algorithm,  and  the  dashed  line  the  recovered  value  using  spatially  guided  reconstruction. 


In  addition  to  improved  qualitative  and  quantitative  accuracy,  the  spatial  prior  algorithm 
reduced  the  reconstruction  time  significantly.  The  full  reconstruction  time  for  the  no-priors 
case,  including  background  optical  property  estimation,  was  just  under  9  minutes  for  both  the 
central  and  superficial  tumor  cases.  These  times  were  reduced  to  less  than  3  minutes  and  90 
seconds  using  spatial  priors.  In  both  cases,  structural  information  guided  the  algorithm  to  a 
convergent  solution  in  far  fewer  iterations  than  the  no-priors  cases. 

4.2  Phantom  results 

Figure  8  shows  fluorescence  yield  images  recovered  from  phantom  data  using  both  no-prior 
and  spatial  prior  image  reconstruction  approaches.  A  qualitative  assessment  of  the  images 
reveals  a  dramatic  benefit  of  the  spatial  prior  on  image  formation.  The  fluorescent  object’s 
borders  are  more  clearly  defined  for  all  concentrations  of  Lutex  when  using  spatial  priors. 
Furthermore,  the  values  of  fluorescence  yield  throughout  the  region  of  interest  are  more 
homogeneous,  and  therefore,  more  similar  to  the  actual  distribution  for  spatially  guided 
reconstructions.  Incorporating  spatial  priors  also  suppresses  edge  artifacts  significantly.  This 
is  most  apparent  in  images  of  phantoms  with  low  Lutex  concentration.  Figure  9  shows  a  full 
scale  image  of  the  0.3125  jliM  phantom.  Artifacts  near  the  boundary  virtually  dominate  the 
no-priors  image  and  represent  the  largest  fluorescent  yield  values  in  that  imaging  domain. 
None  of  these  artifacts  exist  in  the  spatially  guided  image  generated  from  the  same 
fluorescence  emission  data.  The  spatially  guided  image  also  includes  an  easily  discernable 
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fluorescent  object  at  the  correct  location  which  is  difficult  to  define  in  the  no-priors  case. 
These  results  indicate  that  for  low  fluorophore  concentrations  in  this  imaging  domain,  the 
correct  fluorescent  object  would  be  identified  only  if  prior  structural  information  were 
incorporated  in  the  reconstruction. 


Lutex 

Concentration 

5  pM 

EEB 

2.5  pM 

LI  1 

1.25  pM 

LEI 

0.625  pM 

LEI 

0.3125  pM 

LEI 

OpM 

LEI 

Fig.  8.  Recovered  images  of  fluorescence  yield  are  shown  for  varying  concentrations  of  Lutex. 
The  14  mm  diameter  fluorescent  inclusion  was  embedded  in  a  55  mm  diameter  solid  epoxy 
tissue  simulating  phantom.  Images  were  generated  from  the  same  data  using  algorithms  based 
on  no- priors  and  spatial  soft  prior  implementations. 
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Fig.  9.  A  narrower  colorbar-scale  version  of  the  0.3125  pM  Lutex  phantom  images  shown  in 
Fig.  8  further  illustrates  the  improvement  in  image  accuracy  for  the  spatially  guided  algorithm. 

Both  techniques  suffer  from  a  bias  in  the  recovered  fluorescence  yield  which  results  in  a 
positive  value  in  the  region  of  interest  for  the  case  with  no  fluorophore.  This  is  most  likely 
caused  by  bleed-through  of  the  excitation  light  in  the  emission  measurements  coupled  with  an 
imperfect  spectral  fitting  technique  which  results  in  positive  values  for  fluorescence  emission 
intensity  even  in  phantoms  containing  no  Lutex  fluorophore.  Incorporating  additional 
filtering,  such  as  band-pass  source  filtering,  and  optimizing  the  spectral  fitting  routine  should 
address  the  bias  signal  and  decrease  the  potential  for  incorrect  quantification. 

In  addition  to  improving  image  quality  and  fluorescence  yield  quantification,  spatial  prior- 
based  reconstructions  converged  in  substantially  less  time  than  the  no-priors  algorithms. 
Reconstruction  times  ranged  from  50  -  90  seconds  when  not  using  spatial  priors  and 
approximately  22  seconds  for  spatially  guided  reconstructions,  marking  an  improvement  of 
just  over  75%  in  some  cases.  These  numbers  represent  only  the  fluorescence  yield 
reconstructions  and  not  the  recovery  of  background  optical  properties  since  jua  and  jlis’  were 
assumed  known  from  prior  measurements.  In  cases  requiring  the  recovery  of  background 
optical  properties,  the  improvement  in  reconstruction  time  will  be  similar  to  that  previously 
quoted  in  the  simulation  results. 

5.  Discussion 

This  study  introduced  an  effective  method  to  incorporate  MR-derived  tissue  morphology  for 
imaging  fluorescence  yield  at  depth  in  tissue.  The  conceptual  assertion  is  that  diffuse 
tomography  will  be  more  successful  as  an  imaging  modality  when  combined  with  pre-existing 
imaging  systems  which  have  higher  spatial  resolution,  such  as  MRI.  Simulation  and  phantom 
studies  were  used  here  to  validate  the  method  prior  to  implementation  as  a  full  scale  system. 
Using  Lutex  phantom  data,  the  soft  spatial  prior  implementation  improved  qualitative  and 
quantitative  imaging  response  of  fluorescence  yield.  Prior  information  also  suppressed  image 
artifacts  and  more  accurately  represented  the  internal  distribution  of  fluorophore.  In  the  no¬ 
priors  case,  edge  artifacts  dominate  the  image  for  lower  concentrations  of  fluorophore,  and 
provide  a  misleading  interpretation  of  the  internal  distribution  of  the  fluorescent  agent.  The 
phantom  studies  represented  simple  distributions  of  optical  properties  and  fluorophore 
concentration  and  improvements  in  image  quality  were  still  substantial.  It  is  expected  that 
soft  prior  implementations  will  benefit  imaging  performance  to  an  even  greater  extent  in 
complex  tissue  domains,  an  assertion  born  out  in  the  simulation  studies. 

A  simulated  breast  mesh  derived  directly  from  an  MR  image  of  a  human  breast  served  as 
the  test  bed  for  a  complex  tissue  domain.  The  simulations  demonstrated  no  ability  to  recover 
the  internal  distribution  of  the  complicated  domain  without  the  use  of  spatially  guided 
reconstructions.  Qualitatively,  fluorescence  yield  images  generated  without  spatial  priors  had 
little  resemblance  to  the  target  domain  and  completely  disregarded  the  18mm  simulated 
cancer  region,  as  evidenced  by  1-D  cross-sectional  plots  of  fluorescence  yield.  Breakdown  of 
the  images  in  the  no-priors  case  is  likely  due  to  poorly  recovered  background  optical 
properties  as  well  as  the  complexity  of  the  fluorescence  yield  distribution  itself.  Certainly, 
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part  of  the  improvement  in  fluorescence  yield  image  accuracy  can  be  attributed  to  improved 
images  of  background  optical  properties. 

In  these  studies,  it  was  assumed  that  the  fluorophore  distribution  correlates  directly  to  the 
fatty,  fibro-glandular,  and  tumor  tissue  layers,  which  themselves  exhibit  no  significant  intra¬ 
region  heterogeneity,  in  a  given  imaging  domain.  While  this  assumption  appears  reasonable 
based  upon  the  biology  of  the  tissue  for  endogenous  chromophores  [30,  31],  further  studies  of 
the  uptake  of  Lutex  in-vivo  are  required  to  determine  a  correlation  between  tissue  type  and 
fluorophore  localization.  Should  this  assumption  prove  to  be  unrealistic,  the  “soft”  spatial 
prior  approach  offers  some  latitude  in  terms  of  correctly  identifying  structural  prior 
information.  The  algorithmic  implementation  groups  tissue  regions  together  in  a  region- 
specific  regularization  and  allows  individual  nodes  in  those  regions  to  update  independently. 
As  opposed  to  “hard”  prior  approaches  where  nodal  values  in  a  given  region  are  assumed 
homogenous,  a  soft  prior  technique  may  recover  positive  fluorescent  objects  not  directly 
encoded  in  the  spatial  prior  information.  This  is  a  subject  of  further  investigation. 

The  experimental  system  introduced  here  couples  directly  into  a  Philips  3T  MRI  to 
provide  simultaneous  MR  and  NIR  fluorescence  imaging.  Simultaneous  imaging  simplifies 
co-registration  of  the  MR  image  with  the  NIR  reconstruction  domains  and  reduces  overall 
acquisition  time.  Optimization  of  this  system  for  in  vivo  imaging  is  underway  for  both  small 
animal  and  human  breast  imaging. 

6.  Conclusion 

A  spatially  guided  image  reconstruction  implementation  based  on  prior  knowledge  of  tissue 
morphology  was  shown  to  provide  significant  improvements  in  fluorescence  yield  recovery  in 
complicated  tissue  volumes  and  to  be  highly  beneficial  for  simple  domains.  Specifically,  both 
phantom  and  simulation  results  demonstrated  dramatic  improvements  in  recovery  and 
quantification  of  features  in  the  fluorescence  distribution.  Structural  guidance  also  reduced 
image  reconstruction  time  substantially.  A  newly  developed  experimental  system  couples 
full-volume  deep-tissue  fluorescence  spectroscopy  capabilities  directly  into  the  MR  bore  for 
simultaneous  MR-NIR  fluorescence  image  acquisition.  The  results  presented  here  show 
promise  for  this  approach  in  all  cases  and  tissue  volumes  considered.  Extensions  of  this  study 
are  underway  to  determine  the  effect  of  incorrectly  identifying  the  structural  prior,  especially 
in  cases  where  the  MR  images  produce  false  negative  or  false  positive  readings.  A  full 
characterization  of  the  imaging  limits  of  the  experimental  system  will  further  complement  the 
results  presented  here.  It  is  clear  that  incorporating  anatomical  features  derived  from  MR 
images  in  DOFT  image  reconstruction  will  improve  sensitivity  to  lower  concentrations  of 
fluorophore,  qualitative  accuracy,  and  fluorescence  yield  quantification  in-vivo. 
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A  multichannel  spectrally  resolved  optical  tomography  system  to  image  molecular  targets  in  small 
animals  from  within  a  clinical  MRI  is  described.  Long  source/detector  fibers  operate  in  contact 
mode  and  couple  light  from  the  tissue  surface  in  the  magnet  bore  to  16  spectrometers,  each 
containing  two  optical  gratings  optimized  for  the  near  infrared  wavelength  range.  High  sensitivity, 
cooled  charge  coupled  devices  connected  to  each  spectrograph  provide  detection  of  the  spectrally 
resolved  signal,  with  exposure  times  that  are  automated  for  acquisition  at  each  fiber.  The  design 
allows  spectral  fitting  of  the  remission  light,  thereby  separating  the  fluorescence  signal  from  the 
nonspecific  background,  which  improves  the  accuracy  and  sensitivity  when  imaging  low 
fluorophore  concentrations.  Images  of  fluorescence  yield  are  recovered  using  a  nonlinear 
reconstruction  approach  based  on  the  diffusion  approximation  of  photon  propagation  in  tissue.  The 
tissue  morphology  derived  from  the  MR  images  serves  as  an  imaging  template  to  guide  the  optical 
reconstruction  algorithm.  Sensitivity  studies  show  that  recovered  values  of  indocyanine  green 
fluorescence  yield  are  linear  to  concentrations  of  1  nM  in  a  70  mm  diameter  homogeneous  phantom, 
and  detection  is  feasible  to  near  10  pM.  Phantom  data  also  demonstrate  imaging  capabilities  of 
imperfect  fluorophore  uptake  in  tissue  volumes  of  clinically  relevant  sizes.  A  unique  rodent  MR  coil 
provides  optical  fiber  access  for  simultaneous  optical  and  MR  data  acquisition  of  small  animals.  A 
pilot  murine  study  using  an  orthotopic  glioma  tumor  model  demonstrates  optical-MRI  imaging  of  an 
epidermal  growth  factor  receptor  targeted  fluorescent  probe  in  vivo.  ©  2008  American  Institute  of 
Physics.  [DOI:  10.1063/1.2919131] 


I.  INTRODUCTION 

Imaging  fluorescent  molecular  targets  to  characterize  tis¬ 
sue  pathology  in  vivo  is  an  important  objective  with  broad 
implications  for  drug  development  in  small  animals  and 
clinical  diagnosis.  Given  the  availability  of  targeted  molecu¬ 
lar  imaging  agents  for  animal  research,  several  small  animal 
fluorescence  tomography  scanners  have  been  developed  to 
provide  volumetric  images  of  fluorescence  activity.1-6  Other 
studies  have  focused  on  the  systematic  development  and  as¬ 
sessment  of  diffuse  optical  fluorescence  tomographic 
(DOFT)  imaging  techniques  in  tissue  volumes  relevant  to 
human  imaging,  if  targeted  probes  earn  clinical  approval.7-11 
Most  previous  work  involving  human  breast  imaging  has 
been  completed  using  tissue  simulating  phantoms,  though 
fluorescence  tomography  images  of  human  breast  using  the 
nontargeted  fluorophore  indocyanine  green  (ICG)  have  very 
recently  been  published.7 

Pursuing  noninvasive  optical  molecular  imaging  at  depth 
implies  that  the  interrogating  photons  are  measured  only  at 
the  tissue  surface,  making  the  image  inversion  problem  se- 
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verely  underdetermined  and  ill  posed,  even  with  dense 
source  and  detector  configurations.  The  highly  scattered  pho¬ 
ton  held  limits  the  resolution  of  diffuse  optical  tomographic 
(DOT)  images  and  quantitative  interpretation  can  be  diffi¬ 
cult.  However,  researchers  have  found  that  using  tomogra¬ 
phic  data  acquired  at  discrete  near  infrared  (NIR)  wave¬ 
lengths  and  incorporating  the  spectral  features  of  dominant 
tissue  chromophores  in  the  imaging  algorithm  improves  im- 
aging  performance.  In  general,  greater  quantitative  accu¬ 
racy  is  expected  as  the  number  of  wavelengths  increases; 
however,  image  resolution  is  still  relatively  poor  compared  to 
many  clinical  imaging  systems.  Additional  information  may 
be  used  to  guide  image  recovery  by  noting  that  variation  in 
tissue  optical  properties  can  be  loosely  associated  with  the 
tissue’s  anatomical  structure,  providing  an  opportunity  to  in¬ 
corporate  anatomical  information  from  highly  resolved  con¬ 
ventional  imaging  systems  in  the  recovery  of  images  of  op¬ 
tically  relevant  molecules.  In  practice,  standard  magnetic 
resonance  (MR)  images  have  been  used  to  guide  diffusion- 
based  image  reconstruction  of  endogenous  tissue 
chromophores  and  recently  of  fluorescent  molecules. 
An  imaging  system  that  combines  fully  resolved  spectra  of 
transmitted  or  emitted  light  and  simultaneously  acquired 
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FIG.  1.  (Color  online)  Diagram  of  the  MRI-coupled  16  channel  spectros¬ 
copy  system  for  fluorescence  tomography.  Long,  bifurcated  spectroscopy 
fibers  couple  contact  mode  source-detection  directly  into  the  MR  bore  (3  of 
16  fiber  bundles  are  depicted). 

MRI  data  can  provide  a  particularly  rich  data  set  for  a  variety 
of  imaging  applications.  A  design  goal  of  the  system  pre¬ 
sented  here  was  the  development  of  a  multiwavelength  ca¬ 
pable  detection  system  with  a  large  spectral  bandwidth,  pro¬ 
viding  a  platform  to  explore  the  limits  of  broadband 
spectroscopic  tomography  inside  MRI  scanners. 

This  work  describes  a  unique  multichannel 
spectrometer-based  detection  system  for  imaging  fluores¬ 
cence  yield  in  small  animals  and  breast- sized  domains.  The 
system  operates  in  continuous  wave  (cw)  mode,  so  fully 
quantitative  imaging  of  tissues  requires  the  use  of  frequency 
domain  measurements  obtained  with  a  separate  detector 
system.  The  broadband  cw  detection  system  facilitates  the 
use  of  unique  spectral  fitting  preprocessing  techniques  to  de¬ 
couple  the  fluorescence  emission  from  background  contami¬ 
nation.  The  spectroscopic  detection  system  directly  couples 
into  a  Philips  3T  magnet  for  simultaneous  MRI  and  optical 
data  acquisition.  In  this  configuration,  the  highly  resolved 
MR  images  are  used  as  templates  for  imaging  fluorescence 
yield  from  emission  measurements  of  an  injected  fluoro- 
phore.  A  custom  designed  and  manufactured  rodent  coil  pro¬ 
duced  by  Philips  Research,  Hamburg  integrates  the  fiber  op¬ 
tic  array  into  a  small-diameter  radio  frequency  (rf)  pickup 
coil  for  imaging  small  animals  a  3T  MRI. 

II.  SYSTEM  DESIGN 

The  parallel  spectrometer-based  tomographic  imaging 
system,  depicted  in  Fig.  1,  couples  into  a  Philips  3T  MRI 
magnet  and  was  developed  to  be  flexible  enough  to  acquire 
transmission  and  emission  spectra  in  the  NIR.  Major  system 
components  described  below  include  the  spectrograph,  de¬ 
tection  array,  fiber  optic  transmission  of  light  to  and  from  the 
tissue  surface,  the  acquisition  light  source,  and  MRI  coils 
which  incorporate  the  fiber  optic  patient  or  rodent  interfaces. 
Photographs  of  the  system  are  presented  in  Fig.  2. 

A.  Optical  detection  system 

The  optical  detection  system  is  composed  of  16 
Princeton/Acton  Insight:400F  Integrated  Spectroscopy  Sys¬ 
tems  (Acton,  MA)  residing  in  two  custom  designed  wheeled 
carts  (8020,  Columbia  City,  IN).  The  Insight  400F  consists  of 
a  0.3  m  F3.9  imaging  spectrograph  and  a  low  noise,  front 
illuminated  charge  coupled  device  (CCD)  (Pixis  400F) 


FIG.  2.  (Color  online)  The  spectroscopy  system  (a)  is  built  into  two  carts 
which  can  be  wheeled  into  the  MRI  control  room  shown  in  (b).  Four  of  the 
insight  spectrographs  can  be  seen  in  the  cart  on  the  left  in  (a).  Thirteen  meter 
long  fiber  bundles  extend  through  ports  in  the  wall  directly  into  the  MRI 
bore  (c).  The  small  animal  coil  is  shown  in  (d)  with  eight  fiber  bundles 
circumscribing  a  murine  subject. 

cooled  to  -70  °C.  The  1340X400  pixel  CCD  is  vertically 
binned  to  maximize  detector  area/wavelength  providing  a 
binned  pixel  area  of  0.16  mm2.  Manufacturer  specifications 
indicate  a  dark  current  of  0.0025  electrons /pixel  s  and  quan¬ 
tum  efficiencies  of  0.45  at  750  nm  and  0.20  at  950  nm.  Each 
spectrograph  contains  a  motorized  grating  turret  holding  300 
and  1200  1/mm  gratings,  which  when  coupled  to  the  CCD, 
provide  spectral  ranges  of  60  and  300  nm  for  a  single  grating 
position,  respectively.  Both  gratings  are  blazed  at  750  nm  for 
maximum  efficiency  in  the  NIR. 

The  current  system  is  designed  for  cw  imaging  only,  as 
the  CCD  spectrometer  detectors  cannot  measure  rapid  sig¬ 
nals  in  the  hundreds  of  megahertz  frequency  range.  Thus, 
truly  quantitative  imaging  requires  acquiring  data  with  both  a 
frequency  domain  (FD)  system  and  the  CCD-based  spectros¬ 
copy  system.  This  procedure  is  demonstrated  here,  and  the 
future  hardware  integration  of  the  two  is  outlined  in  the 
Discussion. 

B.  Light  collection  and  delivery 

Sixteen  custom  designed  bifurcated  fiber  bundles 
(Zlight,  Latvia)  channel  light  to  and  from  the  imaging  do¬ 
main.  Each  fiber  bundle  is  composed  of  eight  13  m  long  and 
400  jmm  diameter  silica  fibers  which  contact  the  tissue  sur¬ 
face.  Seven  fibers  from  each  bundle  connect  to  the  input 
system  of  each  of  the  spectrographs,  while  the  eighth,  known 
as  a  source  fiber,  branches  off  to  the  source  coupling  system. 
Since  the  detection  branch  of  each  fiber  bundle  is  fixed  to  a 
given  spectrograph,  the  detected  light  path  experiences  no 
fiber-to-fiber  coupling  between  the  tissue  surface  and  spec¬ 
trometer,  minimizing  losses  and  providing  fully  parallel  ac¬ 
quisition.  Details  of  the  light  paths  through  the  spectrograph 
system  are  presented  elsewhere.22  The  light  source  is  sequen- 
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filter  wheel 


FIG.  3.  (Color  online)  Diagram  (top)  and  photograph  (bottom)  of  the  cus¬ 
tom  designed  entrance  optics.  Lenses  collimate  the  incoming  light  for  filter¬ 
ing  with  selectable  interference  or  ND  filters  and  focus  the  light  onto  the 
spectrometer. 

tially  coupled  into  one  of  the  16  source  fibers  using  a  preci¬ 
sion  rotating  stage  (Velmex,  Bloomfield,  NY).  The  detector 
associated  with  the  active  source  fiber  is  deactivated,  result¬ 
ing  in  15  full- spectrum  measurements  acquired  in  parallel  for 
each  of  the  16  source  positions,  though  the  number  of 
source-detector  pairs  may  be  reduced  for  smaller  domains. 

The  seven  fibers  in  the  detection  branch  of  each  fiber 
bundle  are  arranged  in  a  line.  These  bundles  couple  to  the 
spectrographs  via  custom  designed  input  optics  mounted  on 
each  spectrograph,  as  depicted  in  Fig.  3.  An  FI. 4,  25  mm 
focal  length  digital  camera  lens  (Megapixel,  Edmund  Optics, 
Barrington,  NJ)  collects  and  collimates  almost  the  full  cone 
of  light  emitted  from  the  detector  fibers  [numerical  aperture 
(NA)  =  0.37].  The  collimated  light  passes  through  a  fully  au¬ 
tomated  six-position  filter  wheel  containing  two  long  pass 
interference  filters  with  650  and  720  nm  cut-ons  (Omega, 
Brattleboro,  VT)  for  fluorescence  emission  acquisition.  Ad¬ 
ditionally,  two  neutral  density  (ND)  filters  with  optical  den¬ 
sities  (ODs)  of  1  and  2  reside  in  each  filter  wheel  to  increase 
the  dynamic  range  of  transmission  imaging.  Filter  wheel  po¬ 
sitions  are  automatically  adjusted  during  image  acquisition 
depending  on  the  type  of  acquisition.  A  25  mm  diameter, 
60  mm  focal  length  NIR  achromat  (Thorlabs,  Newton,  NJ) 
focuses  the  collimated  and  filtered  light  onto  the  input  slit. 
The  F#  matching  of  fiber  to  spectrograph  utilizes  the  high 
NA  of  the  fibers  to  minimize  the  number  of  fibers  in  the 
bundle.  The  spectrograph  entrance  slits  are  fully  opened  to 
maximize  photon  collection,  and  the  magnified  image  of  the 


fiber  array  at  the  slit  plane  defines  an  effective  slit  width  of 
just  over  1  mm.  Thus,  the  system  operates  at  resolutions  of 
about  2.2  and  11.2  nm  using  the  1200  and  300  1/mm  grat¬ 
ings,  respectively.  The  spectral  resolution  can  be  increased  at 
the  cost  of  throughput  by  reducing  the  input  slit  width. 

Usable  dynamic  range  of  the  16  bit  CCD  chip  is  over  2.8 
orders  of  magnitude  assuming  a  minimum  signal  of  100 
counts  above  background.  Controlling  the  camera  exposure 
times  between  0.01  and  120  s  adds  an  additional  4  orders  of 
magnitude  for  a  total  dynamic  range  approaching  7  orders  of 
magnitude.  This  range  applies  to  fluorescence  emission  mea¬ 
surements  which  are  prefiltered  using  long  pass  filters  in  the 
single- wheel  automated  filter  selector  described  above.  Since 
transmission  mode  measurements  are  not  prefiltering  with 
long  pass  filters,  the  ND  filters  may  be  used  selectively  to 
further  enhance  the  detection  performance,  resulting  in  a  to¬ 
tal  effective  dynamic  range  of  9  orders  of  magnitude.  Maxi¬ 
mum  camera  exposure  times  are  limited  by  the  desired  total 
acquisition  time  and  may  be  increased  to  further  extend  the 
dynamic  range  in  cases  for  which  total  acquisition  time  is  not 
limited. 

C.  Light  sources 

The  imaging  method  for  a  given  acquisition  determines 
which  light  source  to  incorporate  into  the  imaging  sequence. 
Currently  available  sources  in  use  include  a  high  power  tung¬ 
sten  white  light  source  for  broadband  transmission  measure¬ 
ments  and  a  bank  of  laser  diodes.  Since  the  focus  of  the  work 
presented  here  is  fluorescence  imaging,  data  for  the  pre¬ 
sented  examples  were  acquired  using  a  690  nm  cw  laser  di¬ 
ode  (Applied  Optronics,  South  Plainfield,  NJ)  to  excite  the 
fluorophore. 

D.  Patient/animal  interface 

A  unique  MRI  rodent  coil  designed  in  collaboration  with 
Philips  Research  Europe  (Hamburg,  Germany)  features  16 
access  holes  and  nylon  set  screws  to  accommodate  the  spec¬ 
troscopy  fibers  in  a  circular  array.  While  small  magnet  bore 
animal  systems  are  restricted  to  birdcage  coil  configurations, 
larger  bore  commercial  clinical  MR  scanners  allow  the  use 
of  solenoid  coil  designs.  The  Brfield  direction  of  this  coil 
must  be  oriented  orthogonal  to  the  main  B0  field  axis.  For 
small  animal  imaging  at  3T,  the  solenoid  coil  provides  high 
sensitivity  and  good  field  of  view  coverage.  Calculations 
of  electromagnetic  field  and  coil  characteristics  were  per¬ 
formed  from  simulations  using  a  commercially  available  EM 
program  (FEKO),  which  is  based  on  the  method  of  moments 
(Fig.  4). 

The  solenoid  coil  has  an  open  inner  diameter  of  70  mm 
and  is  built  from  8  mm  wide  strip  conductors  wound  around 
a  fiber  glass  cylinder.  The  coil  support  cylinder  is  made  of 
glass  fiber  and  is  mechanically  fixed  to  a  sidewall.  The  par¬ 
allel  windings  are  connected  as  shown  in  Fig.  5  and  the  gap 
between  individual  strip  conductors  is  8  mm,  providing 
ample  space  to  accommodate  the  spectroscopy  fibers.  Cir¬ 
cumscribing  the  fiber  glass  coil  support  cylinder  is  a  remov¬ 
able  polycarbonate  cylinder  containing  nylon  set  screws  to 
affix  the  spectroscopy  fibers. 
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FIG.  4.  (Color  online)  (a)  Relative  transversal  magnetic  field  of  the  solenoid 
coil.  The  isocenter  is  set  to  0  dB  and  the  arrow  shows  the  B0  direction,  (b) 
Sensitivity  profiles  along  the  coil  axis  (perpendicular  to  the  main  field  B0) 
for  different  distances  from  the  isocenter. 

The  overall  inductance  of  the  solenoid  is  Ls- 1024  /zH  at 
127  MHz.  Sixteen  equidistant  splits  with  nonmagnetic  ca¬ 
pacitors  (ATC  100B)  are  introduced  along  the  conductor  to 
avoid  current  inhomogeneities  (propagation  effects).  During 
the  transmit  phase,  the  coil  is  detuned  by  three  independent 
parallel  detuning  circuits  distributed  along  the  solenoid, 
which  make  the  coil  transparent  for  the  B1  transmit  field. 
Optimized  blocking  radio  frequency  chokes  prevent  rf  leak¬ 
age  and  provide  a  high  Q  factor  measured  at  QL=600  for  the 
unloaded  coil.  A  low  noise  preamplifier  is  directly  connected 
to  the  solenoid,  and  optimal  noise  matching  is  performed  via 
a  low  loss  pi  network. 

For  mouse  brain  imaging,  a  custom  designed  molded 
mouse  bed  is  used,  also  shown  in  Fig.  5.  The  bed  contains 
fiber  optic  access  holes  to  provide  stable  positioning  of  eight 
spectroscopy  fibers  around  the  mouse  head,  reducing  the  op¬ 
tical  data  set  to  a  maximum  of  56  source-detector  pairs. 

Larger  volumes  and  patient  breast  imaging  capabilities 
are  provided  by  a  fiber  optic  ring  attached  to  a  commercial 
3T  MRI  breast  coil  (MRI  Devices,  Waukesha,  WI),  depicted 
in  Fig.  6.  The  current  design  requires  manual  fiber  position¬ 
ing  using  set  screws,  however,  spiral  ring  and  parallel 
plate  geometries  may  be  implemented  for  more  reliable 
positioning. 

E.  Automated  acquisition 

The  system  is  operated  using  a  Dell  desktop  personal 
computer  running  Windows  XP  Professional.  Image  acquisi¬ 
tion  is  automated  using  custom  programs  written  in  LAB  VIEW 
(National  Instruments,  Austin,  TX)  developed  with  the  SI- 


(a)  (b) 

FIG.  5.  (Color  online)  (a)  Photograph  of  the  coil  layout  showing  connected 
strip  conductors  and  holes  to  provide  access  for  optical  fibers.  The  custom 
designed  mouse  bed  accommodates  the  spectroscopy  fiber  bundles  and  fits 
into  the  rodent  coil  for  simultaneous  MR  and  optical  acquisition  (b). 


FIG.  6.  (Color  online)  The  patient  interface  is  a  circular  array  of  fibers  that 
couples  into  a  standard  MR  breast  coil  (a)  although  parallel  plate  geometries 
are  also  under  consideration  (b). 

Toolkit  camera/spectrometer  drivers  produced  by  Rcubed, 
LLC  (Lawrenceville,  New  Jersey).  USB  cables  connect  all 
camera/spectrometer/filter  wheel  units  to  the  computer,  and 
the  motorized  source-coupling  stage  is  connected  via  a  serial 
cable.  A  screenshot  of  the  primary  acquisition  control  pro¬ 
gram  is  provided  in  Fig.  7.  The  user  interface  is  designed  to 
minimize  technician  training.  Imaging  modes  include  fluo¬ 
rescence,  transmission,  raw  data  only,  system  calibration, 
and  basis  spectra  acquisition,  and  once  a  mode  is  selected, 
only  relevant  options  pertaining  to  that  imaging  mode  appear 
on  the  screen.  Exposure  times  may  be  manually  set  by  the 
user  or  automatically  determined  using  an  optimization  rou¬ 
tine  which  performs  test  exposures  to  calculate  ideal  expo¬ 
sure  times  for  image  acquisition.  Additional  options  include 
real-time  spectrum  calibration  and  spectral  fitting  and  auto¬ 
matic  neutral  density  filtering  to  preferentially  decrease  filter 
OD  as  a  function  of  source-detector  distance  for  measure¬ 
ments  not  involving  fluorescence  emission. 

III.  FLUORESCENCE  TOMOGRAPHY  IMAGE 
FORMATION 

The  MRI-guided  finite  element  method  (FEM)  image  re¬ 
construction  implementation  is  described  extensively  in  a 
previous  publication.8  In  general,  the  approach  used  here  is 
as  follows. 

(1)  Acquire  FD  data  [photomultiplier  tube  (PMT)  detection 
system]  at  the  excitation  wavelength  and  use  these  data 
to  reconstruct  for  optical  properties,  /Jiax  and  /jlJ  . 

(2)  Acquire  FD  data  (PMT  detection  system)  at  the  emis¬ 
sion  wavelength  and  use  these  data  to  reconstruct  for 
optical  properties,  /zam  and  /ismr . 

(3)  Acquire  cw  transmission  data  by  using  the  fluorescence 
excitation  laser  source  and  the  spectrometer  detection 
system. 

(4)  Acquire  cw  fluorescence  emission  data,  calibrate  these 
data  using  the  data  acquired  in  part  (3),  and  use  the 
reconstructed  optical  properties  to  recover  fluorescence 
yield,  a  product  of  the  fluorophore’s  quantum  efficiency 
7]  and  its  absorption  coefficient,  jmaf(r). 

Currently,  frequency  domain  data  required  for  parts  (1) 
and  (2)  are  collected  using  a  separate  DOT  imaging  system 
used  in  clinical  exams  and  described  in  Ref.  23.  PMT-based 
frequency  domain  capabilities  are  being  integrated  into  the 
spectrometer  system  carts  to  provide  frequency  domain  and 
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FIG.  7.  (Color  online)  Front  panel  of 
the  spectroscopy  system’s  image  ac¬ 
quisition  program. 


cw  spectroscopic  detection  in  a  single  system.  These  modi¬ 
fications  are  outlined  in  the  Discussion  section. 

MR  images  simultaneously  acquired  with  optical  data 
are  used  to  generate  the  finite  element  domain  for  optical 
tomography  image  recovery.  MR  images  are  segmented  us¬ 
ing  thresholding,  region-growing,  and  manual  image  ma¬ 
nipulation  in  materialize  mimics  software  (Ann  Arbor,  MI). 
Source/detector  fiber  optic  positions  are  registered  with  ref¬ 
erence  to  MR  sensitive  fiducials  attached  to  the  patient  or 
animal  fiber  interface.  Exported  mask  files  are  used  to  gen¬ 
erate  FEM  meshes  compatible  with  the  optical  tomography 
modeling  and  image  reconstruction  software  described  in 
Ref.  8 

Two  methods  are  used  to  incorporate  tissue  structural 
information  derived  from  simultaneously  acquired  MR  im¬ 
ages  into  the  optical  tomography  reconstruction  algorithm. 
One  approach  structures  the  inverse  regularization  matrix 
based  on  segmented  MR  images.  This  is  known  as  a  “soft” 
priors  approach  since  each  node  independently  updates,  al¬ 
lowing  the  recovery  of  tumor  regions  not  explicitly  seg¬ 
mented  from  the  MR  image.  Alternatively,  the  “hard  priors” 
approach  homogenizes  each  segmented  region  and  is  there¬ 
fore  unyielding  in  its  application  of  spatial  guidance.24,25 

In  the  soft  priors  approach,  spatial  prior  information  is 
incorporated  by  assuming  a  “generalized  Tikhonov”  penalty 
term  that  results  in  a  Laplacian-type  regularization  matrix.  To 
recover  images  of  optical  parameters,  the  difference  between 
measured  fluence  <Emeas  at  the  tissue  surface  and  calculated 
data  <EC  is  minimized.  The  iterative  update  equation  for  the 
optical  properties  /z=(/za,Ac),  where  pia  is  the  absorption  co¬ 
efficient,  k  is  the  diffusion  coefficient,  k=  l/3(jma+ jul/),  and 
juls '  is  the  reduced  scattering  coefficient,  pertains  to  steps  (1) 
and  (2)  above  and  is  given  by 

A  u  =T  JtJ+B  LrLl_1/r(T>meas  -Oc  )  (1) 

LXrLxjn  J  ^  ^  J  J  ™x,trans_m  ^x,trans_ra/ 

for  the  excitation  and  emission  wavelengths  x  and  m.  Once 
these  values  are  know,  the  fluorescence  yield  update  can  be 
calculated, 


A  wv  =  VTJ  +  /3fLTL]  “ 1 7r(<3>  ™eas  -  <) ,  (2) 

where  in  both  update  equations,  J  is  the  Jacobian  matrix 
describing  the  sensitivity  of  boundary  data  to  the  parameter 
of  interest,  and  (3  is  a  fixed  fraction  multiplied  by  the  maxi¬ 
mum  value  on  the  diagonal  of  JTJ.  Spatial  prior  information 
is  introduced  in  the  dimensionless  “filter”  matrix  L  generated 
using  MR  images  segmented  into  appropriate  tissue  regions 
based  on  MRI  contrast. 

The  hard  priors  approach  applies  stricter  constraints  than 
the  regularization-based  implementation  previously  de¬ 
scribed.  Nodal  values  are  locked  together  by  implicitly  as¬ 
suming  that  segmented  regions  in  the  imaging  domain  are 
homogeneous.  In  practice,  the  Jacobian  matrix  is  calculated 
on  a  fully  resolved  mesh  at  each  iteration,  but  then  collapsed 
into  segmented  regions  for  the  inversion  process.  This  dras¬ 
tically  reduces  the  update  parameter  space  from  the  product 
of  number  of  nodes  and  number  of  unknown  properties  to 
the  product  of  number  of  regions  and  number  of  unknown 
properties.  Confidence  in  the  segmentation  is  critical  since 
regions  are  offered  no  spatial  latitude.  Indeed,  the  hard  priors 
approach  is  strictly  a  characterization-based  implementation 
incapable  of  detection  or  uncovering  false  negatives,  and 
therefore  does  not  present  as  an  image  recovery  problem 
per  se. 

IV.  OPTICAL  DATA  PROCESSING 

The  following  outlines  the  data  processing  procedure  for 
all  optical  data  collected  on  the  spectroscopy  system  (A-C) 
as  well  as  additional  calibration  procedures  for  fluorescence 
emission  measurements  (D). 

A.  Baseline/dark  current  correction 

Each  CCD  displays  a  baseline  offset  for  zero  second 
acquisitions  which  must  be  subtracted  from  the  raw  spec¬ 
trum.  Multiple  repetitions  of  the  baseline  offset  spectrum 
were  measured  for  each  spectrometer.  The  median  of  the 
repetitions  is  used  as  the  baseline  offset  spectrum  to  be  sub- 
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traded,  pixel  by  pixel  (assuming  vertical  binning  of  the  CCD 
chip),  from  each  acquired  spectrum.  This  is  done  for  each 
spectrometer. 

A  similar  method  was  used  to  correct  for  the  dark  current 
of  the  CCDs.  Dark  room  acquisitions  were  recorded  for  a 
range  of  exposure  times.  Measured  counts  are  proportional  to 
exposure  time.  The  counts  per  exposure  time  slope  was  de¬ 
termined  per  binned  pixel  for  each  spectrometer  and  is  used 
to  subtract  dark  current  from  recorded  spectra.  After  baseline 
and  dark  current  correction,  spectra  are  converted  to  counts 
per  second. 

B.  Detector  calibration 

A  first  order  correction  is  applied  to  account  for  hetero¬ 
geneity  in  throughput  and  wavelength  dependence  for  all  op¬ 
tical  components  between  the  tissue  and  CCD  detector  (i.e., 
detection  fibers,  input  optics,  spectrometer  optics,  and  CCD 
response).  Detectors  were  arranged  to  circumscribe  a  6  cm 
diameter  cylindrical  Teflon  phantom  with  an  SMA  connector 
attached  to  the  radial  center  of  one  end.  Light  from  a  high 
power  tungsten  white  light  source  was  focused  into  a  fiber 
connected  to  the  centrally  located  SMA  connector.  Spectra 
were  recorded  with  all  spectrometers  (usually  ten  repeti¬ 
tions),  and  interdetector  calibration  factors  were  calculated 
for  every  vertical  pixel  bin  (total  of  1340)  for  each  spectrom¬ 
eter.  This  was  done  for  each  grating  and  grating  position  to 
be  used  during  image  acquisitions.  Calibration  factors  are 
stored  and  used  to  scale  detected  signal  in  a  wavelength  de¬ 
pendent  manner.  This  helps  reduce  the  influence  of  inhomo¬ 
geneities  in  the  CCD  array  and  partially  accounts  for 
throughput  variability  between  detector  channels. 

A  similar  approach  is  used  to  correct  for  OD  filtering. 
Since  the  ND  filters  have  a  wavelength  dependent  response, 
OD  values  for  every  filter  were  experimentally  determined 
using  the  spectroscopy  system.  For  a  given  spectrometer, 
these  values  were  calculated  for  each  CCD  pixel  (vertically 
binned)  for  each  grating/wavelength  range  selection.  In  this 
manner,  a  calibration  file  for  a  given  filter  setting  at  a  given 
grating  and  wavelength  setting  contains  1340  X  16  values  of 
optical  density,  accounting  for  number  of  pixels  by  number 
of  spectrometers. 

C.  Source  calibration 

A  first  order  interfiber  source  strength  calibration  was 
recorded  by  positioning  a  detector  at  the  radial  center  of  the 
Teflon  cylinder  described  above.  The  detected  signal  from  all 
16  source  fiber  positions  was  used  to  determine  source  scal¬ 
ing  factors  for  future  acquisition. 

D.  Fluorescence  emission  spectral  fitting 
and  data-model  calibration 

Fluorescence  emission  light  is  initially  decoupled  from 
the  often  much  more  intense  excitation  signal  using  long 
pass  interference  filters  which  provide  between  5  and  7  OD 
filtering  efficiency  depending  on  the  filter  set  and  the  wave¬ 
length  selectivity  of  the  spectrograph  grating.  Even  with 
these  components  in  place,  light  reaching  the  CCD  retains 
residual  signal  not  associated  with  fluorescence  emission 


.  Measured  Spec. 

-  Fit  Spec. 


-  Excitation  basis 

-  Fluorescence  basis 


FIG.  8.  (Color  online)  A  spectral  fitting  routine  is  used  to  decouple  fluores¬ 
cence  spectra  from  the  measured  signal.  The  spectral  fitting  algorithm  de¬ 
termines  the  relative  contribution  of  the  basis  spectra  [right  column,  (b)  and 
(d)]  which  best  fits  the  measured  spectrum.  Two  examples  are  presented 
above,  both  for  a  70  mm  diameter  homogeneous  phantom  containing  in¬ 
tralipid,  ink,  and  ICG  [500  nM  for  graphs  (a)  and  (b)  and  100  pM  for  graphs 
(c)  and  (d)]. 


from  the  fluorophore  target,  which  can  be  a  significant  com¬ 
ponent  of  the  total  signal,  especially  for  cases  in  which  the 
fluorophore  has  a  low  quantum  yield  or  is  at  low  concentra¬ 
tions  in  the  tissue.  In  order  to  further  decouple  the  true  fluo¬ 
rescence  signals,  previously  recorded  “basis  spectra”  of  the 
residual  nonfluorophore  originating  signal  and  the  pure  fluo¬ 
rescence  signal  are  fitted  to  the  data.  Each  spectrograph  is 
used  to  record  its  own  fluorescence  emission  basis  spectra 
and  residual  excitation  basis  spectra  in  a  controlled  manner. 
All  spectra  are  normalized  to  the  maximum  value.  Fitting  is 
accomplished  by  a  linear  least  squares  algorithm  that  mini¬ 
mizes  the  summation 

N 

S  =  2{j,-[«W  +  W?(V)]}2  (3) 

1=1 

with  respect  to  a  and  b ,  where  yt  is  the  measured  intensity  at 
a  given  wavelength  pixel,  F  and  G  are  the  residual  excitation 
and  fluorescence  basis  spectra,  a  and  b  are  the  coefficients 
recovered  in  the  minimization  procedure,  and  A  is  the  num¬ 
ber  of  wavelength  pixels  per  spectrum.  The  algorithm  deter¬ 
mines  the  amount  of  fluorescence  emission  and  excitation 
cross-talk  in  the  measured  spectrum,  the  sum  of  which  best 
fits  the  data  in  a  least  squares  sense.  This  procedure  is  ap¬ 
plied  to  each  recorded  fluorescence  spectrum,  a  total  of  240 
spectra  for  a  given  acquisition.  The  spectral  fitting  procedure 
is  demonstrated  for  two  different  ICG  concentrations  in  a 
homogeneous  phantom  in  Fig.  8.  In  some  cases,  the  fluores¬ 
cence  peak  has  been  observed  to  shift  to  longer  wavelengths 
in  large  phantom  volumes.  To  account  for  this,  the  minimi¬ 
zation  may  incorporate  the  position  of  the  peak  as  a  free 
parameter.  The  resulting  fluorescence  emission  peak  is  inte- 
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grated  to  provide  a  single  fluorescence  intensity  measure¬ 
ment  for  a  given  source-detector  pair. 

Integrated  fluorescence  emission  measurements  are  cali¬ 
brated  in  the  following  manner: 


^calib 


vmeas 

fh 


^model 

xi 

(F)meas_spec  ’ 
xi 


(4) 


where  the  index  i  indicates  a  single  data  point  or  source- 
detector  pair.  O™odel  is  the  intensity  boundary  data  calculated 
from  the  images  of  jnax  and  /jlsx  recovered  using  the  fre¬ 
quency  domain  measurements  from  the  PMT  clinical  system. 
This  calculated,  or  “model,”  boundary  data  must  correspond 
to  the  spectrometer  fiber  positions,  which  do  not  necessarily 
have  to  be  identical  to  those  used  in  the  frequency  domain 
system.  The  quotient  in  Eq.  (4)  essentially  provides  a  scaling 
factor  for  the  fluorescence  measurements  which  accurately 
scales  the  fluorescence  data  to  the  FEM  model  and  accounts 
for  inter-detection  channel  throughput  and  fiber  coupling  dis¬ 
crepancies.  An  initial  estimate  for  the  iterative  reconstruction 
algorithm  is  determined  using  a  homogeneous  fitting  routine. 
This  procedure  uses  the  bisection  method  to  minimize  the 
data-model  misfit,  assuming  a  homogeneous  distribution  of 
fluorescence  yield. 


V.  SYSTEM  PERFORMANCE 
A.  Repeatability 

An  8.6  cm  diameter  homogeneous  phantom  composed 
of  silicone,  titanium  dioxide,  and  India  ink  was  used  to  mea¬ 
sure  the  repeatability  of  transmission  mode  measurements 
using  the  spectroscopy  system.  Optical  properties  of  the 
phantom  were  approximately  0.004  and  1.91  mm-1  for  the 
absorption  and  reduced  scattering  coefficients,  respectively. 
Measurements  were  repeated  at  each  source-detector  position 
using  the  690  cw  laser  source.  For  a  single  vertically  binned 
pixel  array,  the  average  and  maximum  standard  errors  at  the 
laser  peak  are  1.4%  and  1.8%,  respectively,  if  fiber  positions 
were  not  changed  in  between  measurements.  These  increase 
to  18.3%  and  19.5%  with  fiber  repositioning,  indicating  that 
fiber  coupling  is  the  most  significant  source  of  error  in  these 
measurements.  If  pixels  are  binned  throughout  the  laser 
peak,  the  average  and  maximum  standard  errors  change  to 
0.28%  and  0.37%  for  fibers  remaining  in  contact  with  the 
phantom  and  14.8%  and  16.3%  for  repositioned  fibers. 

Determining  measurement  repeatability  in  fluorescence 
mode  is  less  straightforward  given  the  signal  dependence  on 
fluorophore  concentration,  absorption  spectrum,  and  quan¬ 
tum  yield.  For  this  study,  a  70  mm  diameter  liquid  phantom 
containing  DPBS,  1%  intralipid,  India  ink,  and  10  nM  ICG 
was  used  to  determine  measurement  repeatability.  This  mea¬ 
sure  was  calculated  in  two  ways,  one  considered  only  the 
raw  data  for  a  given  binned  pixel  array  and  resulted  in  a 
mean  standard  error  of  0.7%  and  maximum  standard  error  of 
1.4%  for  the  pixel  array  at  the  fluorescence  peak.  The  second 
measure  was  determined  based  on  the  integrated  intensity 
from  the  full  calibration  and  spectral  fitting  routine,  resulting 
in  an  average  standard  error  of  0.6%  and  maximum  standard 
error  of  1.8%.  It  is  clear  that  fiber  positioning  variability 
would  dominate  the  data  error  for  fluorescence  measure¬ 


ments  in  this  case.  However,  fluorescence  measurement  data 
offer  a  unique  opportunity  to  account  for  fiber  coupling  vari¬ 
ability.  This  is  accomplished  by  calibrating  the  fluorescence 
measurements  to  the  transmission  measurements  in  the  same 
geometry,  as  described  in  Sec.  IV  D.  This  provides  inherent 
stability  to  systematic  error,  though  is  not  unique  to 
spectrometer-based  detection  which  in  and  of  itself  does  not 
necessarily  provide  a  signal  to  noise  ratio  (SNR)  advantage 
over  more  traditional  detection  schemes  used  for  fluores¬ 
cence  tomography.  The  benefit  of  spectrometer  based  detec¬ 
tion  is  improved  separation  of  the  fluorescence  signal  from 
background  contamination,  providing  a  more  accurate  ratio, 
especially  for  low  concentrations  of  fluorophore. 

B.  Sensitivity 

A  70  mm  diameter  liquid  phantom  containing  dulbecco’s 
phosphate  buffered  saline  (DPBS),  1%  intralipid,  and  India 
ink  was  used  to  investigate  the  overall  sensitivity  of  the  sys¬ 
tem  to  ICG  fluorescence.  ICG  dye  dissolved  in  de-ionized 
(DI)  water  was  added  to  the  solution  to  obtain  solutions 
ranging  from  10  pM  to  1  jmM  ICG.  The  optical  properties  of 
the  intralipid/ink  solution  fxa  and  juls'  were  approximately 
0.006  and  1.6  mm-1,  respectively.  The  optical  properties  at 
the  excitation  and  emission  wavelengths  were  recovered  for 
each  concentration  of  ICG  from  data  collected  using  the 
clinical  frequency  domain  system.  Since  the  domain  was 
known  to  be  homogeneous,  these  properties  were  determined 
in  a  homogeneous  fitting  procedure  only. 

Fluorescence  emission  and  excitation  transmission  mea¬ 
surements  were  recorded  for  each  phantom,  with  a  maximum 
allowed  camera  integration  time  of  120  s  applied  to  the  fluo¬ 
rescence  measurements.  The  typical  data  measured  across 
the  emission  spectrum  are  shown  in  Fig.  8(a)  for  a  strong 
fluorescence  signal  and  Fig.  8(c)  for  a  weak  signal.  The  spec¬ 
tral  fitting  procedure  discussed  in  Sec.  IV  D  above  was  used 
to  recover  the  true  fluorescence  signal  and  the  nonspecific 
background  contributions,  as  shown  in  (b)  and  (d),  respec¬ 
tively.  Additionally,  to  quantitatively  compare  the  spectral 
fitting  procedure  to  more  conventional  means  of  fluorescence 
filtering,  the  measured,  un-fitted  spectra  were  integrated  to 
simulate  720  nm  long-pass  filtering.  In  both  cases,  integrated 
values  representing  fluorescence  emission  intensities  were 
calibrated  as  per  Eq.  (4)  and  used  to  determine  homogeneous 
values  of  fluorescence  yield. 

Values  of  fluorescence  yield  recovered  using  the  two 
data  preprocessing  techniques  are  plotted  as  a  function  of 
known  ICG  concentration  in  Fig.  9.  The  linear  fit  shown  in 
the  figure  was  computed  using  the  spectrally  fit  data  with  the 
y-intercept  forced  to  zero  and  indicates  a  strong  linear  corre¬ 
lation  of  recovered  fluorescence  yield  and  ICG  concentration 
(R2=0.99).  At  concentrations  above  1  nM,  the  fluorescence 
yield  values  calculated  using  data  that  was  filtered  only,  with 
no  spectral  fitting,  very  closely  match  the  spectrally  fit  re¬ 
sults.  Fluorescence  signals  produced  at  these  fluorophore 
concentrations  dominate  the  detected  signal,  a  finding  con¬ 
sistent  with  the  data  shown  in  Fig.  8.  However,  fluorescence 
yield  values  recovered  with  spectrally  fit  data  maintain  the 
linear  relationship  at  lower  concentrations  than  those  recov¬ 
ered  using  data  without  the  spectral  separation  of  back- 
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FIG.  9.  (Color  online)  Values  of  fluorescence  yield  as  a  function  of  ICG 
concentration  in  a  70  mm  liquid  phantom  recovered  using  two  methods  to 
process  the  recorded  data.  One  method  uses  the  spectral-fitting  technique  to 
decouple  background  contamination  across  the  fluorescence  emission  spec¬ 
trum,  while  the  other  simply  integrates  the  measured  spectrum,  as  would  be 
the  case  for  long  pass  filtering  alone.  Values  were  determined  using  a  ho¬ 
mogeneous  fitting  algorithm  for  the  background  optical  properties  as  well  as 
the  fluorescence  yield. 

ground  contamination.  Clearly,  the  two  techniques  diverge  at 
1  nM.  At  this  concentration,  the  residuals  between  the  recov¬ 
ered  value  of  fluorescence  yield  and  the  linear  approximation 
are  0.4%  and  48%  for  the  spectrally  fit  and  filtering-only 
approaches,  respectively.  As  the  fluorophore  concentration 
drops  below  1  nM,  both  techniques  lose  the  consistent  linear 
response  observed  at  higher  concentrations,  though  the 
filtered-only  data  show  even  less  sensitivity  to  changes  in 
fluorophore  concentration.  At  500  pM,  the  residuals  increase 
to  57%  and  192%  for  the  spectrally  fit  and  filtered-only  data 
processing  responses,  respectively.  Below  this  level,  accurate 
quantification  of  fluorescence  activity  is  unlikely  in  this 
phantom  configuration,  however,  the  spectrally  fit  data  still 
show  a  stronger  response  to  changes  in  fluorophore  concen¬ 
tration  down  to  10  pM.  The  residuals  calculated  at  the  lowest 
concentration  measured  were  almost  five  times  larger  with¬ 
out  spectral  fitting,  clearly  indicating  a  more  sensitive,  if  not 
accurate,  response  of  the  spectral  preprocessing  technique.  It 
should  be  noted  that  since  fluorophore  quantum  yield  is  not 
explicitly  known  in  this  solution,  the  calculated  slope  of  the 
linear  regression  does  not  provide  information  on  the  actual 
relationship  between  true  and  recovered  concentrations, 
however,  the  linearity  itself  is  a  critical  measure  of  system 
performance. 


FIG.  10.  (Color  online)  Images  of  fluorescence  yield  for  a  70  mm  diameter 
liquid  phantom  containing  an  ICG  heterogeneity.  The  target-to-background 
concentration  of  ICG  was  3.3:1.  Images  reconstructed  using  spatial  soft 
priors,  depicted  in  (b),  are  qualitatively  more  accurate  and  in  terms  of  re¬ 
covered  contrast  than  images  recovered  without  the  spatial  priors  implemen¬ 
tation  (a). 

1  /x M,  providing  a  total  contrast  of  just  over  3.3:1.  Back¬ 
ground  optical  properties  were  determined  by  imaging  the 
phantom  in  a  separate  frequency  domain  clinical  system  and 
values  at  785  nm  were  used  for  this  experiment.  The  spec¬ 
troscopy  system  was  then  used  to  acquire  excitation  intensity 
and  fluorescence  emission  measurements.  Integrated  intensi¬ 
ties  were  calibrated  to  the  model  as  described  earlier,  and 
images  were  reconstructed  using  an  algorithm  without  spatial 
guidance  as  well  as  the  soft  spatial  priors  implementation. 
The  resulting  images  are  shown  in  Fig.  10.  Tumor- to- 
background  contrast  can  be  deciphered  in  both  images;  how¬ 
ever,  the  spatially  guided  implementation  provides  a  more 
accurate  representation  of  the  imaged  domain.  Inclusion  bor¬ 
ders  are  more  clearly  defined  and  the  image  represents  sub¬ 
stantially  higher  tumor-to-background  contrast. 

D.  Small  animal  imaging 

An  animal  pilot  study  demonstrates  the  ability  of  the 
system  to  image  an  epidermal  growth  factor  receptor 
(EGFR)  targeted  fluorescent  probe  in  mouse  brain  tumors. 
U-251  human  glioma  tumors  were  implanted  intracranially 
in  male  nude  mice  19  days  prior  to  the  imaging  study.  Si¬ 
multaneous  gadolinium  enhanced  T1  MRI  and  fluorescence 
tomography  acquisition  of  the  population  were  completed 
prior  to  IRDye®  800CW  EGF  Optical  Probe  (LI-COR  Bio¬ 
sciences,  Lincoln,  NE)  IV  administration  (1  nM)  and  at  24  h 
intervals  after  injection  over  the  course  of  72  h.  Eight  fibers 
were  used  in  this  study,  and  fiber- tissue  contact  positions 
were  determined  in  reference  to  MRI  sensitive  fiducials  on 
the  molded  mouse  bed.  For  each  imaging  session,  a  single 
coronal  slice  was  segmented  into  regions  using  MATERIALIZE 
MIMICS  software.  The  example  in  Fig.  11  shows  an  MR  coro- 


C.  Phantom  imaging 

A  heterogeneous  liquid  phantom  was  used  to  demon¬ 
strate  imaging  large  volumes  with  imperfect  tumor-to- 
background  uptake.  The  phantom  was  composed  of  DPBS, 
1%  intralipid,  and  India  ink,  resulting  in  background  optical 
properties  of  /xa  =  0.005  mm-1  and  jm/  =  1.4  mm-1.  ICG  dis¬ 
solved  in  DI  water  was  added  to  the  phantom  volume  to 
obtain  a  300  nM  ICG  solution.  A  thin-walled  plastic  cylinder 
was  positioned  between  the  edge  and  center  of  the  phantom 
to  simulate  a  2  cm  diameter  tumor  region.  The  inclusion  con¬ 
sisted  of  the  same  solution  found  in  the  phantom  back¬ 
ground,  although  the  ICG  concentration  was  elevated  to 


FIG.  11.  (Color  online)  A  gadolinium  enhanced  T1 -weighted  image  shows  a 
coronal  slice  of  a  tumor  bearing  mouse  head  (a).  The  image  is  segmented 
into  three  regions  using  mimics  software  (b),  the  largest  region  (violet)  iden¬ 
tifies  the  tissue  outside  of  the  brain,  the  second  largest  region  (gray)  is  the 
brain  itself,  and  the  smallest  region  (light  blue)  marks  the  tumor.  The  seg¬ 
mented  image  is  used  to  generate  an  FEM  mesh  for  spatially  guided  fluo¬ 
rescence  reconstruction. 
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FIG.  12.  (Color  online)  Reconstructed  images  of  fluorescence  yield  for  a 
tumor  bearing  mouse  using  only  the  outer  boundary  for  spatial  priors  (a)  and 
using  spatial  hard  priors  (b). 

nal  slice  of  one  mouse  48  h  after  IRDye®  800CW  EGF  Op¬ 
tical  Probe  administration  and  the  resulting  MR  image.  The 
regionized  mask  was  exported  as  a  bitmap  file  from  which 
a  two  dimensional  mesh  was  generated  for  fluorescence 
reconstruction. 

The  spectral  fitting  technique  described  in  Sec.  IV  D  was 
used  to  generate  integrated  intensity  values  of  fluorescence 
emission  from  basis  spectra  of  EGF  Optical  Probe  fluores¬ 
cence.  Reconstructed  images  with  and  without  the  use  of 
hard  spatial  priors  are  presented  in  Fig.  12.  The  conventional 
diffuse  tomography  approach  makes  no  use  of  the  internal 
tissue  structure  although  the  outer  boundary  of  the  domain 
was  used  in  this  case.  Fluorescence  yield  values  recovered 
using  this  approach  are  highly  surface  weighted,  showing  an 
elevated  region  of  fluorescence  activity  near  the  surface  of 
the  mouse  head,  well  outside  the  actual  location  of  the  tumor. 
The  diffuse,  Gaussian  shape  of  the  elevated  region  is  typical 
of  fluorescent  molecular  tomography  imaging  without  spatial 
guidance. 3,5,26-31  The  use  of  hard  priors,  as  implemented 
here,  transforms  the  imaging  problem  to  one  of  quantifica¬ 
tion  since  the  recovery  of  the  tumor  region  size  and  location 
is  not  the  ultimate  objective.  In  this  case,  the  tumor  region 
outline  is  given  by  a  region-growing  threshold  of  the  en¬ 
hanced  MR  image,  and  the  fluorescence  recovery  algorithm 
estimates  the  fluorescence  uptake  in  that  predefined  region. 
An  important  test  will  be  how  the  hard  priors  implementation 
handles  false  positives  in  healthy  regions,  and  nontumor 
bearing  mice,  introduced  by  incorrect  MR  segmentation,  and 
these  studies  are  ongoing. 

VI.  DISCUSSION 

An  MRI-coupled  optical  tomography  system  has  been 
developed  to  image  fluorescence  yield  in  a  variety  of  tissue 
volumes.  A  unique  rodent  coil  specifically  developed  for  this 
system  provides  high  resolution  MR  images  of  small  ani¬ 
mals.  Through  a  variety  of  spatial  prior  implementations, 
these  images  provide  the  imaging  template  upon  which  the 
fluorescence  activity  is  reconstructed.  The  system’s  high  sen¬ 
sitivity,  low  noise  CCD  detectors  were  shown  to  be  stable 
and  provide  repeatable  measurements,  the  most  significant 
error  originating  from  fiber  coupling  errors.  Fluorescence 
emission  intensity  may  be  scaled  to  the  excitation  laser  trans¬ 
mitted  intensity  to  account  for  these  errors. 

Characterizing  the  system’s  sensitivity  to  fluorophore 
concentration  is  complex,  given  the  dependence  on  concen¬ 
tration,  absorption  spectra,  excitation  filtering,  and  quantum 
yield.  Sensitivity  results  of  a  commonly  used  fluorophore  in 


a  reasonably  sized  phantom  were  presented  and  demon¬ 
strated  a  linear  response  of  the  fluorescence  yield  fit  to  fluo¬ 
rophore  concentration  down  to  1  nM  using  spectral  fitting 
data  preprocessing.  Without  spectral  fitting,  recovered  values 
for  1  nM  and  below  are  inaccurate,  and  the  sensitivity  slope 
begins  to  flatten  at  these  concentrations.  This  reduction  in 
sensitivity  is  also  observed  for  spectrally  fit  data,  although 
the  decrease  in  the  response  slope  is  less  dramatic  and  does 
not  show  up  until  concentrations  of  500  pM.  Below  this 
level,  there  is  still  a  modest  response  to  changing  fluorophore 
concentrations,  although  the  recovered  values  are  clearly  in¬ 
accurate.  This  is  not  caused  by  a  drop  in  measurement  SNR 
per  se ,  as  the  measured  signal  even  at  10  pM  is  strong  and 
stable,  but  may  be  due  to  systematic  bias  in  the  spectral 
fitting  routine  itself,  inaccuracies  in  the  basis  spectra,  or  the 
presence  of  a  contaminating  signal  with  similar  spectral  char¬ 
acteristics  to  ICG. 

Preliminary  imaging  studies  demonstrate  an  ability  to 
image  low  fluorescence  contrast  (3.3:1  ratio)  in  relatively 
large  tissue  volumes  (70  mm  diameter),  contingent  upon  the 
availability  of  spatial  prior  information  provided  by  the  MRI. 
Additionally,  it  was  also  shown  that  quantifying  the  uptake 
of  a  targeted  fluorophore  in  a  mouse  glioma  model  is  feasible 
and  dramatically  benefits  from  the  inclusion  of  structural  in¬ 
formation,  through  the  synergy  of  MR  imaging  and  fluores¬ 
cence  estimation  of  the  regionized  image.  A  larger  mouse 
population  study  is  currently  underway  to  further  character¬ 
ize  in  vivo  imaging  performance,  and  quantify  the  variation 
in  estimated  uptake  of  fluorophores  in  tumors. 

The  work  presented  here  required  the  use  of  a  separate 
frequency  domain  imaging  system  to  obtain  background  op¬ 
tical  properties  and  calibrate  the  fluorescence  data  to  the 
algorithmic  model.  This  is  inconvenient  at  best  and  limits 
many  of  the  advantages  of  the  MRI-coupled  imaging  system. 
Importantly,  the  lack  of  FD  imaging  capabilities  integrated 
into  the  imaging  system,  background  optical  properties  may 
not  be  acquired  without  repositioning  the  subject.  Integrating 
PMT-based  frequency  domain  detection  into  the  spectros¬ 
copy  system  is  underway.  The  rotating  source  coupling  stage 
was  constructed  with  15  PMT  detectors  in  a  design  identical 
to  that  previously  reported.23  In  this  configuration,  the  source 
branch  of  the  spectroscopy  fibers  will  serve  as  both  light 
delivery  and  pickup  channels  for  the  PMT’s. 

The  unique  system  introduced  here  provides  rich  spec¬ 
tral  information  in  a  variety  of  diffuse  optical  imaging  modes 
including  broadband  NIR  transmission  and  fluorescence  at 
present.  In  fluorescence  mode,  the  spectrograph  system  of¬ 
fers  several  advantages  over  filtered  intensity  measurements 
as  spectrally  resolved  detection  provides  exceptional  wave¬ 
length  selectivity  for  excitation  filtering.  The  excitation  con¬ 
tamination  and  nonspecific  background  can  be  dramatically 
suppressed  to  improve  the  ability  to  quantify  low  level  fluo¬ 
rescence,  as  was  shown  here.  The  system  may  also  be  used 
in  the  future  to  acquire  emission  data  from  multiple  fluoro¬ 
phores  simultaneously,  contingent  upon  the  individual  fluo¬ 
rescence  peaks  being  resolved  well  enough  for  the  spectral 
fitting  technique  to  extract  the  contribution  from  each 
fluorophore. 


Downloaded  16  Jul  2008  to  129.170.203.126.  Redistribution  subject  to  AIP  license  or  copyright;  see  http://rsi.aip.org/rsi/copyright.jsp 


064302-10  Davis  et  al. 


Rev.  Sci.  Instrum.  79,  064302  (2008) 


ACKNOWLEDGMENTS 

This  work  was  funded  by  the  National  Institutes  of 
Health  Grants  No.  ROl  CA109558,  ROl  CA69544,  and  U54 
CA 105480,  Philips  Research  Hamburg,  as  well  as  Depart¬ 
ment  of  Defense  Breast  Cancer  predoctoral  fellowship 
BC051058. 

1N.  Deliolanis,  T.  Lasser,  D.  Hyde,  A.  Soubret,  J.  Ripoll,  and  V. 
Ntziachristos,  Opt.  Lett.  32,  382  (2007). 

2H.  Meyer,  A.  Garofalakis,  G.  Zacharakis,  S.  Psycharakis,  C.  Mamalaki,  D. 
Kioussis,  E.  N.  Economou,  V.  Ntziachristos,  and  J.  Ripoll,  Appl.  Opt.  46, 
3617  (2007). 

3S.  Y.  Patwardhan,  S.  R.  Bloch,  S.  Achilefu,  and  J.  P.  Culver,  Opt.  Express 
13,  2564  (2005). 

4J.  Ripoll,  R.  B.  Schulz,  and  Y.  Ntziachristos,  Phys.  Rev.  Lett.  91,  103901 
(2003). 

5R.  B.  Schulz,  J.  Ripoll,  and  V.  Ntziachristos,  IEEE  Trans.  Med.  Imaging 
23,  492  (2004). 

6G.  Zacharakis,  J.  Ripoll,  R.  Weissleder,  and  V.  Ntziachristos,  IEEE  Trans. 
Med.  Imaging  24,  878  (2005). 

7  A.  Corlu,  R.  Choe,  T.  Durduran,  M.  A.  Rosen,  M.  Schweiger,  S.  R. 
Arridge,  M.  D.  Schnall,  and  A.  G.  Yodh,  Opt.  Express  15,  6696  (2007). 
8S.  C.  Davis,  H.  Dehghani,  J.  Wang,  S.  Jiang,  B.  W.  Pogue,  and  K.  D. 
Paulsen,  Opt.  Express  15,  4066  (2007). 

9 A.  Godavarty,  M.  J.  Eppstein,  C.  Zhang,  and  E.  M.  Sevick-Muraca, 
Radiology  235,  148  (2005). 

10  A.  Godavarty,  A.  B.  Thompson,  R.  Roy,  M.  Gurfmkel,  M.  J.  Eppstein,  C. 
Zhang,  and  E.  M.  Sevick-Muraca,  J.  Biomed.  Opt.  9,  488  (2004). 

11  A.  Godavarty,  C.  Zhang,  M.  J.  Eppstein,  and  E.  M.  Sevick-Muraca,  Med. 
Phys.  31,  183  (2004). 

12 R.  Choe,  A.  Corlu,  K.  Lee,  T.  Durduran,  S.  D.  Konecky,  M.  Grosicka- 
Koptyra,  S.  R.  Arridge,  B.  J.  Czerniecki,  D.  L.  Faker,  B.  Chance,  M.  A. 
Rosen,  and  A.  Yodh,  Med.  Phys.  32,  1128  (2005). 

13  A.  Corlu,  R.  Choe,  T.  Durduran,  K.  Lee,  M.  Schweiger,  S.  R.  Arridge,  E. 

M.  Hillman,  and  A.  Yodh,  Appl.  Opt.  44,  2082  (2005). 

14 A.  Corlu,  T.  Durduran,  R.  Choe,  M.  Schweiger,  E.  M.  Hillman,  S.  R. 
Arridge,  and  A.  G.  Yodh,  Opt.  Lett.  28,  2339  (2003). 


15  A.  Li,  Q.  Zhang,  J.  Culver,  E.  Miller,  and  D.  Boas,  Opt.  Lett.  29,  256 
(2004). 

16  S.  Srinivasan,  B.  W.  Pogue,  B.  Brooksby,  S.  Jiang,  H.  Dehghani,  C. 
Kogel,  W.  A.  Wells,  S.  P.  Poplack,  and  K.  D.  Paulsen,  Technol.  Cancer 
Res.  Treat.  4,  513  (2005). 

17 S.  Srinivasan,  B.  W.  Pogue,  S.  Jiang,  H.  Dehghani,  C.  Kogel,  S.  Soho,  J. 
J.  Gibson,  T.  D.  Tosteson,  S.  P.  Poplack,  and  K.  D.  Paulsen,  Acad.  Radiol. 
13,  195  (2006). 

18 S.  Srinivasan,  B.  W.  Pogue,  S.  Jiang,  H.  Dehghani,  and  K.  D.  Paulsen, 
Appl.  Opt.  44,  1858  (2005). 

19  B.  Brooksby,  S.  Jiang,  C.  Kogel,  M.  Doyley,  H.  Dehghani,  J.  B.  Weaver, 
S.  P.  Poplack,  B.  W.  Pogue,  and  K.  D.  Paulsen,  Rev.  Sci.  Instrum.  75,  5262 
(2004). 

20M.  Guven,  B.  Yazici,  X.  Intes,  and  B.  Chance,  Phys.  Med.  Biol.  50,  2837 
(2005). 

21 X.  Intes,  C.  Maloux,  M.  Guven,  B.  Yazici,  and  B.  Chance,  Phys.  Med. 
Biol.  49,  N155  (2004). 

22 B.  W.  Pogue,  Z.  Li,  C.  Carpenter,  A.  Laughney,  V.  Krishnaswamy,  S.  C. 
Davis,  S.  Jiang,  and  K.  D.  Paulsen,  International  Society  for  Optical 
Engineering  (SPIE)  BiOS  at  Photonics  West,  San  Jose,  California,  2008 
(unpublished). 

23 T.  O.  McBride,  B.  W.  Pogue,  S.  Jiang,  U.  L.  Osterberg,  and  K.  D.  Paulsen, 
Rev.  Sci.  Instrum.  72,  1817  (2001). 

24 H.  Dehghani,  B.  W.  Pogue,  S.  Jiang,  B.  Brooksby,  and  K.  D.  Paulsen, 
Appl.  Opt.  42,  3117  (2003). 

25 M.  Schweiger  and  S.  R.  Arridge,  Phys.  Med.  Biol.  44,  2703  (1999). 

26  V.  Ntziachristos,  C.  H.  Tung,  C.  Bremer,  and  R.  Weissleder,  Nat.  Med.  8, 
757  (2002). 

27  E.  Graves,  J.  Ripoll,  R.  Weissleder,  and  V.  Ntziachristos,  Med.  Phys.  30, 
901  (2003). 

28  V.  Ntziachristos,  E.  A.  Schellenberger,  J.  Ripoll,  D.  Yessayan,  E.  Graves, 
A.  Bogdanov,  Jr.,  L.  Josephson,  and  R.  Weissleder,  Proc.  Natl.  Acad.  Sci. 
U.S.A.  101,  12294  (2004). 

29  E.  E.  Graves,  R.  Weissleder,  and  V.  Ntziachristos,  Current  Molecular 
Medicine  4,  419  (2004). 

30 G.  Zacharakis,  J.  Ripoll,  R.  Weissleder,  and  Y.  Ntziachristos,  IEEE  Trans. 
Med.  Imaging  24,  878  (2005). 

31 E.  E.  Graves,  D.  Yessayan,  G.  Turner,  R.  Weissleder,  and  V.  Ntziachristos, 
J.  Biomed.  Opt.  10,  0440191  (2005). 


Downloaded  16  Jul  2008  to  129.170.203.126.  Redistribution  subject  to  AIP  license  or  copyright;  see  http://rsi.aip.org/rsi/copyright.jsp 


Spectral  Distortion  in  Diffuse  Molecular  Luminescence 


Tomography  in  Turbid  Media 

Scott  C.  Davis1,  Brian  W.  Pogue1,  Stephen  B.  Tuttle1,  Hamid  Dehghani1,2,  Keith  D. 
Paulsen1 

1.  Thayer  School  of  Engineering,  Dartmouth  College,  Hanover,  New  Hampshire  03755 

2.  School  of  Physics,  University  of  Exeter,  Exeter,  UKEX4  4QL 
Scott.C.Davis@dartmouth.edu,  Brian.  W.Pogue@dartmouth.edu 

ABTRACT 

The  influence  of  tissue  optical  properties  on  the  shape  of  near-infrared  (NIR) 
fluorescence  emission  spectra  propagating  through  multiple  centimeters  of  tissue-like 
media  is  investigated.  Experimentally  measured  fluorescence  emission  spectra  measured 
in  6  cm  homogeneous  tissue-simulating  phantoms  shows  dramatic  spectral  distortion 
which  results  in  emission  peak  shifts  of  up  to  60  nm  in  wavelength.  Measured  spectral 
shapes  are  highly  dependent  on  the  photon  path-length  and  the  highly  scattered  photon 
field  in  the  NIR  amplifies  the  wavelength-dependent  absorption  of  the  fluorescence 
spectra.  Simulations  of  the  peak  propagation  using  diffusion  modeling  describe  the 
experimental  observations  and  confirm  the  path-length  dependence  of  fluorescence 
emission.  Spectral  changes  are  largest  for  longer  path-length  measurements,  and  this  will 
be  most  dominant  in  human  tomography  studies  in  the  NIR.  Spectrally  resolved  or  multi¬ 
wavelength  band-pass  measurements  are  required  to  detect  these  changes,  and  may  be 
essential  to  interpret  such  effects,  which  would  otherwise  be  attributed  to  erroneous 
intensity  measurement.  This  phenomenon  is  analogous  to  beam  hardening  in  x-ray 
tomography,  which  can  lead  to  image  artifacts  without  appropriate  compensation.  The 


peak  shift  toward  longer  wavelengths,  and  therefore  lower  energy  photons,  observed  for 
NIR  luminescent  signals  propagating  through  tissue  may  readily  be  described  as  a  beam 
softening  phenomenon. 


INTRODUCTION 


Diffuse  spectroscopy  and  imaging  of  tissue  using  luminescent  signals  from 
molecular  reporters  has  become  a  major  area  of  interest  in  recent  years1'4.  Much  of  this 
research  is  focused  on  exploiting  novel  biochemical  probes,  developing  new  detection 

r  o 

strategies,  improving  image  recovery  methods  '  ,  and  coupling  luminescence  detectors  to 
existing  imaging  systems9,10.  Image  recovery  or  quantitative  spectroscopy  of  tissue 
requires  accurate  light  transport  modeling  which  is  now  a  fairly  mature  field  of  research 
which  has  seen  substantial  growth  since  the  early  1990s11.  While  fluorescence  molecular 
imaging  and  spectroscopy  systems  for  use  in  small  animals4,10,12'16  have  evolved  to  the 

i  n 

point  of  commercial  availability  and  clinical  trials  in  humans  are  underway  ,  some  of  the 
more  subtle  complexities  of  the  signal  acquisition  remain  to  be  examined.  One  important 
consideration  that  has  not  been  investigated  in  detail,  especially  in  the  near-infrared 
(NIR),  is  the  interaction  between  the  remission  spectrum  of  the  reporter  and  the 
intervening  tissue  through  which  the  light  transport  occurs  prior  to  detection.  In 
particular,  changes  in  spectral  emission  can  be  detected  which  are  intimately  linked  to  the 
tissue  optical  properties  and  the  path-length  of  signal  travel  involved.  In  this  study,  a 
systematic  evaluation  of  NIR  spectral  shift  has  been  completed  using  simulations  and 
tissue  phantoms. 

As  interest  in  fluorescence  spectroscopy  for  tissue  diagnosis  grew,  researchers 
began  to  develop  methods  to  compensate  for  spectral  distortion  due  to  tissue  optical 
properties  in  order  to  recover  intrinsic  fluorescence  spectra.  Original  work  reported  by 
Wu  et  al.18,  Durkin  et  al.19  and  Richards-Kortum  et  al.20,  examined  photon  migration, 
Kubelka-Munk  and  exponential  models,  respectively,  to  extract  intrinsic  autofluorescence 


signals  from  the  measured  distorted  spectrum  emitted  from  tissue.  In  1996,  Durkin  et 
al.  compared  several  modeling  approaches  and  determined  that  a  partial  least  squares 
method  yielded  an  accurate  spectral  correction.  Analytical  expressions  derived  by 
Gardner  et  al.  were  used  to  extract  spectra  measured  through  tissue  samples  on  the 
order  of  1  cm.  In  2001,  Muller  et  al.23  provided  a  comprehensive  investigation  of  the 
effects  of  absorption  and  scattering  on  intrinsic  fluorescence  extraction  based  on  a  photon 
migration  model. 

All  of  these  efforts  were  focused  on  tissue  spectroscopy  in  the  visible  spectrum 
where  autofluorescence  is  high  and  the  significant  photon  attenuation  restricts  the 
distances  over  which  light  signals  can  be  measured  to  a  few  centimeters.  The  relatively 
low  attenuation  in  the  near- infrared  (NIR)  allows  measurable  light  penetration  over  10 
cm,  though  elastic  scattering  ensures  most  photons  will  change  direction  within  1  mm  of 
entering  the  tissue.  Photon  propagation  in  the  NIR  is  commonly  modeled  using  the 
diffusion  approximation  to  the  radiative  transport  equation.  While  Patterson  and  Pogue24 
developed  a  general  construct  for  modeling  fluorescence  propagation  through 
homogeneous  tissues  based  on  diffusion  theory,  they  did  not  consider  the  spectral 
distortion  in  the  detected  emissions  directly.  Some  investigators25,26  have  recognized  the 
inherent  depth  information  contained  in  the  spectral  distortion  and  have  demonstrated  the 
ability  to  localize  fluorescent  layers  in  tissue  simulating  phantoms  by  calculating  the  ratio 
of  emitted  intensity  at  different  wavelengths. 

With  the  development  of  NIR  diffuse  tomography  for  imaging  endogenous  tissue 
contrast  through  centimeters  of  tissue,  extensions  of  the  theoretical  framework  to 
incorporate  fluorescence  detection  have  been  sought  to  enhance  the  disease-to-normal 
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tissue  contrast  and  target  molecular  processes  using  fluorescent  contrast  agents  '  . 
Systems  have  been  developed  to  collect  fluorescence  intensities  from  an  array  of  source- 
detector  positions  surrounding  the  tissue  of  interest  and  model-based  reconstruction 
techniques  have  been  used  to  recover  fluorescence  activity  by  matching  the  modeled  and 
measured  data  using  optimization  techniques.  Most  experimental  systems  employ  one  or 
more  band-pass  or  long-pass  filters  to  separate  the  excitation  signal  from  the  fluorescence 
emission,  but  rarely  has  the  capability  to  resolve  the  fluorescence  spectrum  been 
incorporated.  Recently,  it  was  shown10,34  that  a  spectrally  resolved  fluorescence 
tomography  scanner  that  employs  spectrometer  detection  can  record  fluorescence  spectra 
emitted  from  tissue.  Ignoring  the  spectral  changes  in  deep  tissue  imaging  may  have 
implications  for  quantitative  accuracy  of  the  recovered  images,  especially  for  systems 
that  use  relatively  large  emission  wavelength  ranges.  Since  wavelength-dependent 
attenuation  of  the  fluorescence  emission  increasingly  distorts  the  measured  spectrum  as 
photon  path-length  increases,  the  effect  is  expected  to  become  more  pronounced  when 
imaging  through  larger  tissue  domains,  such  as  the  human  breast. 

This  work  investigates  NIR  fluorescence  emission  distortion  using  homogeneous 
turbid  phantoms  and  diffusion  based  modeling.  The  model  platform,  NIRFAST,  has 
been  developed  previously34,35  using  the  Finite  Element  Method  (FEM)  which  allows 
consideration  of  arbitrarily  shaped  domains  with  heterogeneous  optical  properties.  Here, 
the  model  system  is  extended  to  calculate  spectrally  resolved  emission  spectra  at  discrete 
wavelengths  based  on  fluorescence  spectra  of  dilute  solutions  and  wavelength-dependent 
optical  properties  of  tissue.  Experimentally  measured  spectra  are  compared  with  spectra 
generated  using  the  model  system.  Phantom  geometries  are  representative  of  the 


dimensions  to  be  encountered  in  tomographic  imaging  of  fluorescence  activity  in  vivo 
and  the  potential  implications  that  spectral  distortion  has  for  tomographic  imaging  of 
fluorescence  are  discussed. 

EXPERIMENTAL  DETAILS 

Phantom  design 

Phantoms  designed  to  simulate  the  optical  properties  of  commonly  imaged  tissue 
types,  typically  adipose  and  fibro-glandular  breast  tissue,  were  used  to  investigate  the 
spectral  distortion  of  the  detected  emission  of  a  distributed  fluorescent  drug.  Two 
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techniques  were  used  to  produce  tissue-simulating  phantoms  ’  .  Liquid  phantoms  were 
composed  of  water  or  phosphate  buffered  solution  (PBS)  and  intralipid,  which  is  a  fatty 
suspension  that  scatters  light  much  like  real  tissue.  Typically,  a  1%  intralipid  solution 
can  be  used  to  mimic  tissue  scattering  in  the  NIR.  The  fluorescent  contrast  agent  may 
then  be  added  in  the  desired  concentration.  Absorbing  suspensions  such  as  India  ink  or 
whole  blood  may  also  be  added  to  provide  appropriate  amounts  of  optical  absorption 
beyond  that  introduced  by  the  fluorophore,  though  these  were  not  used  in  this  study. 
Though  simple  to  prepare  and  flexible  in  terms  of  composition,  the  container  walls 
introduce  unwanted  heterogeneity  and  light  piping  effects. 

Gelatin  phantoms  composed  of  gelatin  extracted  from  porcine  skin  (Sigma 
Aldrich),  TiCL,  blood  and  the  exogenous  contrast  agent  can  be  made  to  closely  mimic 
heterogeneous  tissue  without  artificial  boundaries.  The  mixing  procedure  used  here 
began  by  stirring  40g  of  Type  A  porcine  skin  gelatin  (Sigma  Aldrich)  into  500ml  water  or 
PBS.  Once  dissolved,  the  solution  was  heated  to  about  40  C  in  a  microwave  oven  until 


the  solution  was  fully  transparent  (approximately  1  minute).  Immediately  after  heating, 
the  solution  was  stirred  for  40  minutes.  In  cases  where  scattering  was  required,  0.4g  of 
Ti02  powder  (Sigma  Aldrich)  was  slowly  added  at  the  beginning  of  the  stirring 
procedure.  It  is  critical  that  the  media  appear  homogeneously  cloudy  with  few  visible 
aggregates  of  Ti02  after  20  minutes  of  stirring.  Towards  the  end  of  the  40  minutes,  the 
fluorescent  drug  was  added  while  the  solution  is  still  viscid.  Once  mixed  well,  the 
solution  was  poured  into  a  mold  lined  with  petroleum  jelly  and  refrigerated  for  at  least 
one  hour.  The  cured  phantom  slid  out  of  the  mold  easily.  Photographs  of  liquid  and 
gelatin  phantoms  used  in  this  study  are  shown  in 
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Fig.  1. 

The  fluorescent  drug  investigated  was  Lutetium  Texaphyrin  (LuTex),  a  water 
soluble  dye  developed  as  a  photodynamic  therapeutic  sensitizer.  It  is  a  texaphyrin  based 
molecule  with  lutetium  in  the  center  to  induce  a  large  triplet  state  splitting38.  The 
absorbance  and  fluorescence  emission  peaks  in  dilute  solution  occur  at  approximately 
735  nm  and  750  nm,  respectively.  Fig.  2  shows  a  normalized  fluorescence  emission 
spectrum  of  LuTex  in  a  dilute  solution  of  DI  (deionized)  water  measured  by  a  standard 
fluorometer  (Yvon  Jobin). 

Experimental  system 

The  measurement  system  was  a  fluorescence  tomography  scanner  designed  to 
detect  low  levels  of  light  emitted  from  thick  tissue  samples  (up  to  10  cm).  A  detailed 
description  of  this  instrument  can  be  found  in  Davis  et  al.10.  The  imaging  array  used  in 
the  studies  reported  here  consisted  of  16  source/detector  fibers  surrounding  a  cylindrical 
tissue  phantom  as  shown  in  Fig.  3.  Each  fiber  is  bifurcated  with  one  branch  terminating 
on  a  rotary  stage  which  serially  couples  light  from  a  690  nm  CW  laser  diode  into  any  one 
of  the  fibers.  The  other  branch  couples  light  from  the  phantom  surface  into  one  of  16 
Insight:400F  Spectroscopy  units  with  cooled  CCD  detection  (Princeton  Instruments, 
Acton  MA).  Automated  filter  wheels  filter  the  detected  light  before  it  enters  the 
spectrograph  slit.  In  the  studies  presented  here,  720  nm  long  pass  filters  were  used  to 
suppress  the  excitation  light  in  the  measured  fluorescence  emission  signal.  In  this  study, 
only  a  single  source  fiber  was  illuminated  and  fluorescence  measurements  recorded  at 


various  locations  around  the  tissue  phantom  and  thus  the  system  was  not  used  to  acquire 
tomographic  data. 


Modeling  Light  Propagation 

Near-infrared  light  propagation  in  tissue  is  often  modeled  by  the  diffusion 
approximation,  a  simplification  of  the  radiative  transport  equation  applicable  in  media 
dominated  by  photon  scattering.  Modeling  fluorescence  excitation  and  emission  in  tissue 
is  accomplished  through  a  system  of  diffusion  equations  which  describe  the  photon  field 
produced  by  incident  excitation  light  and  fluorescence  emission  from  a  distributed 
fluorophore.  Specifically,  the  excitation  source,  q0(r) ,  produces  an  excitation  fluence 

rate,  O^r),  throughout  a  medium  with  optical  properties  pax (r)  and  p'sx(r)  representing 

the  absorption  and  reduced  scattering  coefficients  respectively.  The  field  of  excitation 
photons  drives  the  fluorescence  emission  which  propagates  through  the  tissue  subject  to 
optical  properties  pam  (r)  and  psm  (r)  at  the  emission  wavelengths.  The  coupled 

equations  describing  this  process  are  presented  for  the  continuous  wave  case  which 
matches  the  instrumentation  in 
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Fig.  1  27  32: 

-  V  •  kx  (r,  4  )  VOx  (r,  /lx  )  +  (//ax  (r,  4  ))Ox  (r,  Ax)  =  q0(ryAr)  ( 1 ) 

- V •  ^ (r, 4 )VOyz (r, 4 )  +  (Km (r, A!m ))0/ (r,A!J  =  Ox (r)r/(A!m )<uaf (r)  (2) 
where  subscripts  x  and  m  represent  the  excitation  and  emission  fluence  at  wavelengths  Xx 

and  respectively  and  the  diffusion  coefficient  is  k :xm- - - — T - •  The  source 
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ax,m  A sx,m  ) 

term  in  Eq.  (2)  is  a  product  of  the  excitation  photon  field  and  the  fluorescence  yield, 
defined  as  7(4)Aa/(r)  >  which  itself  is  a  product  of  the  fluorophore’s  quantum  efficiency 
rj  and  its  absorption  coefficient,  juaf  (r) ,  at  the  excitation  wavelength. 

Unlike  the  excitation  field,  which  is  assumed  to  arise  from  a  light  source  at  a 
single  excitation  wavelength,  the  emission  spectrum  of  the  fluorophore  can  cover  several 
hundred  nanometers.  To  account  for  this  spectral  bandwidth  in  the  numerical  model,  Eq. 
(2)  is  discretized  over  a  wavelength  range  indicated  by  the  index  i.  Though  the  excitation 
field  and  fluorophore  absorption  coefficient  in  the  source  term  of  Eq.  (2)  do  not  depend 
on  emission  wavelength,  the  fluorescence  source  strength  varies  with  wavelength  based 
on  the  shape  of  the  emission  spectrum  of  the  fluorophore.  This  is  introduced  into  the 
model  through  a  wavelength-dependent  fluorescence  quantum  yield,  r/(/ lm ) .  The 
fluorescence  emission  propagates  through  the  tissue  subject  to  wavelength-dependent 
optical  properties,  juam(r,lm)  and  /usm(r,Xm),  which  must  be  modeled  to  accurately  describe 
the  fluorescence  spectrum  measured  at  the  tissue  surface.  In  practice,  this  is 


accomplished  by  first  calculating  the  excitation  field  in  Eq.  (1),  and  then  solving  for  the 
emission  field  in  Eq.  (2)  at  wavelength  Xm  for  each  i. 

Type  III  boundary  conditions  (also  known  as  Robin  or  mixed)  are  used  to 
describe  the  fractional  loss  of  photons  at  the  tissue-air  interface.  The  flux  leaving  the 
external  boundary  is  described  by: 


(%)+2An  •  k  (£)V<D  (£)  =  0 


(3) 


where  %  is  a  point  on  the  external  boundary,  and  A  is  derived  from  Fresnel’s  law  and 
depends  upon  the  relative  refractive  index  (RI)  mismatch  between  tissue,  Q,  and  air; 


A  = 


2/(l-R0)-l+  cos0c 
1  -  cos  6C  2 


(4) 


where  0C  is  the  angle  at  which  propagation  from  within  the  domain  undergoes  total 
internal  reflection  at  the  boundary  and 
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The  FEM  is  used  to  discretize  the  domain  of  interest  for  numerical  modeling  of 
the  light  propagation  which  was  performed  with  the  NIRFAST  software  package 
developed  for  forward  and  inverse  diffuse  tomographic  imaging  of  tissue.  A  complete 
description  of  the  numerical  implementation  can  be  found  elsewhere  ’  ’  ’  . 


The  emission  spectrum  of  LuTex  was  modeled  from  700  nm  to  850  nm  in 
intervals  of  at  most  10  nm.  Some  intervals  were  smaller  given  the  availability  of 
information  at  additional  wavelengths.  The  extinction  spectrum  of  LuTex  was  measured 
directly  using  a  Cary  50  UV-Vis  spectrophotometer  (Varian,  Inc.,  Palo  Alto  CA). 
Absorption  coefficients  at  each  wavelength  were  calculated  as  a  sum  of  the  constituent 
chromophores,  usually  oxygenated  and  deoxygenated  hemoglobin,  water,  and  the 
exogenous  agent  of  interest,  however,  since  whole  blood  was  not  used  in  the  phantom 
experiments,  only  water  and  LuTex  were  considered  in  the  numerical  portion  of  the 
study.  Published  values  for  the  extinction  spectrum  of  water,  compiled  by  Prahl40,  were 
used  to  calculate  the  absorbing  contribution  of  water.  Scattering  properties  of  tissue  in 
the  NIR  are  commonly  modeled  using  an  empirical  approximation  integrated  into  the 
NIRFAST  software41. 

RESULTS  AND  DISCUSSION 

Phantom  results 

Emission  spectra  measured  through  homogeneous  gelatin  and  water  phantoms 
containing  5  pM  LuTex  and  no  scattering  media  are  presented  in  Fig.  4.  The  phantoms 
were  60  mm  diameter  cylindrical  shapes  with  a  single  source  and  multiple  detectors 
around  the  outer  surface  in  one  plane.  Data  presented  are  for  different  source-detector 
position  on  the  boundary.  As  the  photon  path-length  increases,  the  influence  of  the  media 
changes  the  emission  spectra  modestly  and  thus  the  emission  peak  is  not  significantly 
distorted  when  measured  through  these  phantoms. 


In  turbid  phantoms  of  the  same  size,  the  emission  spectra  change  dramatically. 
LuTex  emission  spectra  measured  for  a  range  of  source-to-detector  distances  around  a 
cylindrical  turbid  phantom  are  presented  in  Fig.  5.  Data  for  spectra  measured  through 
both  gelatin  and  liquid  phantoms  are  shown  and  illustrate  the  significance  of  the  peak 
distortion  through  several  centimeters  of  turbid  media.  Source-detector  distances  were 
1.1,  3.3,  4.9,  and  5.9  cm  for  the  spectra  included  in  the  figure.  Within  1  cm,  the 
measured  fluorescence  spectrum  is  similar  to  the  dilute  sample,  though  a  secondary  peak 
is  evident  around  800  nm  in  both  phantoms.  As  the  source-detector  distance  increases, 
the  spectra  are  changed  more  dramatically.  The  secondary  peak  becomes  increasingly 
prominent  and  as  the  propagation  distance  increases,  the  entire  spectrum  settles  to  a 
single  peak  at  820  nm. 

Photon  path-length  may  be  altered  by  source-detector  geometry  or  by  changes  in 
scattering  properties.  Fig.  6  demonstrates  the  influence  of  intralipid  concentration  on  the 
measured  emission  spectrum.  The  spectral  changes  observed  as  result  of  increasing 
intralipid  concentration  are  similar  to  those  observed  with  increasing  source-detector 
distance,  as  expected. 

Diffusion  modeling 

When  modeling  photon  propagation  through  the  phantoms,  it  was  assumed  that 
the  dominant  absorbers  were  water  (100%)  and  LuTex.  Extinction  spectra  of  water  and 
LuTex  were  used  to  calculate  the  absorption  coefficients  at  discrete  wavelengths  across 
the  range  covered  by  the  emission  spectrum.  These  values  are  plotted  as  a  function  of 
wavelength  in  Fig.  7  and  illustrate  the  large  change  in  absorption  due  mostly  to  the 


fluorophore’s  absorbance.  These  values  were  used  to  determine  the  fluorescence 
intensity  at  each  wavelength  throughout  the  emission  spectrum. 

Emission  intensities  determined  from  the  diffusion  model  at  discrete  10  nm 
wavelength  intervals  from  700  to  850  nm  are  shown  in  Fig.  8  for  the  same  source- 
detector  positions  measured  through  phantoms  and  presented  in  Fig.  5.  Experimental 
data  from  intralipid  phantoms  are  also  plotted  for  comparison  (blue  line)  and  the  red 
dotted  line  represents  the  fluorescence  peak  in  a  dilute  solution.  All  spectra  are 
normalized  to  their  peaks.  Significantly,  the  general  trends  recorded  in  the  phantom 
experiments  are  also  observed  for  the  numerical  model.  At  the  detector  nearest  the 
source,  about  1  cm,  the  measured  fluorescence  emission  peak  is  shifted  about  10  nm  to 
the  red  and  a  small  secondary  peak  is  visible  at  810  nm,  similar  to  what  was  observed  in 
the  phantom  data.  This  corresponds  to  a  very  small  bump  in  the  absorption  spectrum  at 
800  nm,  which  emerges  as  a  dip  in  the  emission  spectrum  and  creates  the  secondary  peak. 
Increased  source-detector  distances  result  in  more  substantial  changes  to  the  emission 
spectrum,  and  an  amplification  of  the  influence  of  the  small  elevation  in  absorption  at 
800  nm.  The  fluorescence  distortion  recorded  at  the  detector  positioned  just  over  3  cm 
from  the  source  results  in  a  peak  at  790  nm  and  a  more  pronounced  dip  at  800  nm  is 
observed,  producing  a  stronger  secondary  peak.  The  original  peak  at  750  nm  is  almost 
entirely  absorbed.  At  longer  source-detector  distances,  the  emission  peak  settles  at  830 
nm,  similar  to  the  phantom  results,  though  the  increase  of  absorption  at  800  nm  is  more 
pronounced  for  the  simulation  results. 

Qualitative  assessment  of  Fig.  8  reveals  reasonable  agreement  between  model  and 
experimental  observations  in  turbid  phantoms  in  terms  of  emission  spectrum  shape, 


especially  for  detector  locations  far  from  the  source.  Pearson  correlation  coefficients 
between  the  model  and  data  were  calculated  as  0.50,  0.90,  0.97  and  0.94  for  the  source- 
detector  distances  of  1.1,  3.3,  4.9,  and  5.9  cm,  respectively.  The  presence  of  tissue 
scattering  increases  the  path-length  of  photons  propagating  through  the  tissue,  which  in 
turn  amplifies  the  influence  of  the  tissue’s  absorption  spectrum.  The  significance  of  the 
peak  distortion  also  indicates  that  the  primary  source  of  fluorescence  detected  at  the  far- 
source  detectors  is  generated  close  to  the  source  excitation  source,  where  the  excitation 
field  is  most  intense,  and  then  propagates  through  the  tissue.  This  is  an  expected  result 
since  photon  penetration  depth  at  the  excitation  wavelength  is  lower  than  within  the 
emission  wavelength  range  due  to  higher  absorption. 

Implications  for  fluorescence  tomography  in  deep  tissue 

The  extent  of  the  spectral  distortion  observed  will  depend  on  the  fluorophore 
used,  the  tissue  chromophore  composition,  and  the  excitation  and  measurement 
geometry.  Other  drugs  used  for  imaging  may  experience  less  distortion  than  that 
associated  with  LuTex.  However,  the  observed  phenomenon  has  implications  for 
quantitative  fluorescence  tomographic  imaging.  Clearly,  the  results  presented  here 
indicate  that  spectrally  resolved  or  band-pass  filtering  is  preferred  over  long-pass  filtering 
of  fluorescence  signals. 

To  demonstrate  the  extent  to  which  spectrally  unresolved  data  impacts 
fluorescence  imaging  in  deep  tissue,  a  simulated  example  using  an  86mm  circular  test 
domain  to  represent  a  coronal  slice  of  a  human  breast  is  considered.  To  create  a  realistic 
test  domain,  tissue  optical  properties  were  calculated  directly  using  the  extinction  spectra 


of  the  dominant  endogenous  chromphores  and  typical  tissue  concentrations  of  those 
chromophores.  LuTex  was  assumed  to  be  the  exogenous  fluorrophore  and  extinction 
coefficients  of  both  endogenous  and  exogenous  chromophores  used  in  this  example  are 
shown  in  Fig.  9  along  with  the  spatial  distributions  of  the  chromophores  for  the  coronal 
slice  under  consideration.  In  this  case,  two  small  Lutex  heterogeneities  were  included  at 
contrasts  of  just  over  3:1  and  2.5:1  over  the  background.  A  large  heterogeneity  with 
contrast  in  hemoglobin  and  water  was  also  added.  Scattering  properties  were  held 
constant  in  the  domain. 

Noise-free  data  were  generated  using  the  forward  model  described  above  and 
signal  contamination  due  to  excitation  cross-talk  and  tissue  autofluorescence  was 
ignored.  Images  of  fluorescence  yield  were  recovered  using  diffusion-based  optimization 
reconstruction  techniques  described  extensively  elsewhere  34,35.  Identical  reconstruction 
algorithms  were  used  to  consider  data  detected  using  one  of  two  measurement 
approaches.  The  first  approach  assumed  that  measured  data  was  fully  wavelength 
resolved.  The  intensity  recorded  at  the  tissue  surface  at  a  single  wavelength  was 
extracted  and  optical  properties  at  that  wavelength  were  used  in  the  model-based 
reconstruction  algorithm.  The  second  approach  simulated  an  experimental  system  using 
only  long-pass  filters.  Since  spectral  selectivity  within  the  emission  spectrum  is 
impossible  in  this  configuration,  the  measured  spectrum  at  each  detector  was  integrated. 
Tissue  optical  properties  were  chosen  to  match  those  at  750  nm,  the  fluorescence 
emission  peak  in  dilute  solution.  Although  the  spatial  heterogeneity  of  these  optical 
properties  was  assumed  to  be  known  exactly  at  this  wavelength,  the  lack  of  spectral 


resolution  does  not  account  for  the  variation  of  optical  properties  across  the  integrated 
spectrum. 

Reconstructed  images  using  both  techniques  are  shown  alongside  the  target  image 
of  fluorescence  yield  in  Fig.  10.  The  only  difference  between  the  two  approaches  arises 
from  the  handling  of  the  detected  fluorescence  spectra.  Images  reconstructed  using 
spectrally  resolved  data  at  a  single  wavelength  show  accurate  recovery  of  fluorescence 
activity  in  the  domain.  Recovered  values  in  the  contrast  enhanced  regions  are  slightly  off 
target  values,  though  this  is  expected  given  the  reported  trade-off  between  enhanced 
region  size  and  recovered  contrast  .  Overall,  quantification  and  localization  are 
excellent.  Certainly,  this  is  to  be  expected  given  the  simple  geometry  and  noise-free  data. 
However,  the  data-model  misfit  error  introduced  by  integrating  the  full  spectrum  is  too 
much  for  the  imaging  algorithm  and  complete  breakdown  in  imaging  performance  is 
observed,  even  for  this  relatively  simple  domain.  The  data  provide  no  ability  to  localize 
or  quantify  fluorescence  contrast  enhancement.  The  failure  of  the  algorithm  using  these 
data  indicates  that  the  long-pass  filtering  approach  is  intractable  for  experimental 
imaging,  making  spectrally-resolved  detection  imperative  for  accurate  image  recovery. 
The  extent  to  which  this  applies  can  be  investigated  by  varying  the  wavelength  range 
over  which  the  spectrum  is  integrated  to  determine  the  widest  filter  band  width  at  which 
image  artifacts  are  insignificant. 

The  dramatic  influence  of  chromophore  absorption  on  the  fluorescence  spectrum 
leaving  the  tissue  clearly  indicates  that  the  spectral  shape  contains  information  related  to 
the  location  of  the  emitting  source.  It  is  reasonable  to  postulate  that  this  information  can 
be  exploited  to  improve  tomographic  fluorescence  imaging  in  vivo,  presuming  spatial 


distributions  of  tissue  chromophore  concentrations  are  known  a  priori  and  the  measured 
emission  spectrum  contains  insignificant  background  signal.  Similar  approaches  have 
been  used  in  bioluminescence  imaging  of  small  animals43,44  and  may  easily  be  adapted  to 
fluorescence  image  reconstruction. 


SUMMARY 

The  problems  of  spectral  shift  are  well  known  in  many  imaging  applications  such 
as  x-ray  computed  tomography,  where  high  atomic  number  materials  can  cause  a  larger 
than  expected  attenuation  of  the  longer  wavelength  photons,  and  shift  the  spectrum  to 
higher  energies,  a  phenomenon  known  as  beam  hardening45.  Luminescent  peaks 
propagating  through  tissue  in  the  NIR,  on  the  other  hand,  exhibit  shifting  towards  lower 
energy  photons,  a  phenomenon  which  softens  the  photon  field.  The  principles  in  both 
cases  are  quite  similar  in  that  the  spectrum  to  be  detected  is  not  attenuated  equally  as  it 
traverses  the  medium. 

A  dramatic  fluorescence  emission  peak  shift  in  the  NIR  was  demonstrated  in 
tissue  simulating  phantoms,  and  this  effect  is  especially  pronounced  when  the  tissue  path- 
lengths  over  which  the  light  signal  is  measured  become  larger.  Thus,  the  effect  will 
likely  have  a  greater  impact  on  imaging  human  organs  than  in  small  animal  applications. 
The  influence  of  tissue  absorption  was  shown  to  produce  the  spectral  distortion,  an  effect 
amplified  by  increasing  photon  path-lengths  caused  by  tissue  photon  scattering  or 
extended  source-to-detector  distances.  Emission  spectra  escaping  the  tissue  phantom  can 
be  modeled  reasonable  well  using  a  finite  element  formulation  of  diffusion  theory, 
especially  for  source-detector  distances  over  2  cm.  These  spectral  changes  should  be 


considered  for  fluorescence  tomography  imaging  through  several  centimeters  of  tissue. 
Use  of  spectrally  resolved  detection  allows  quantification  of  this  change,  and  may  be  the 
only  reliable  way  to  track  intensity  changes  which  would  otherwise  appear  erroneous. 
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FIGURE  CAPTIONS 


Fig.  1 .  Photographs  of  phantoms  used  in  this  study.  A  liquid  intralipid-based  phantom  (a) 
readily  allows  changing  of  the  fluorophore  concentration  while  gelatin  phantoms 
containing  TiOEUB  scatterer  shown  in  (b)  eliminate  the  need  for  an  external  container. 
Gelatin  phantoms  without  scattering  material  were  also  used  to  assess  the  emission 
spectra  without  scatter  (c). 

Fig.  2.  Normalized  LuTex  fluorescence  emission  as  measured  in  dilute  concentration 
(2pg/ml)  in  deionized  water. 

Fig.  3.  A  diagram  of  the  experimental  system  with  16  spectrometers  (Princeton 
Instruments  Inc,  Acton  MA)  is  provided  in  (a)  with  a  photograph  shown  in  (b).  Each 
spectrograph  is  coupled  to  a  fiber  optic  which  surrounds  a  tissue  phantom  in  a  circular 
geometry  (c).  A  homogeneous  teflon  calibration  phantom  is  pictured  here  with  the  16 
fibers  coupled  to  the  exterior  in  a  circle. 

Fig.  4.  Fluorescence  emission  of  5uM  LuTex  measured  in  the  multi-spectrometer  system 
in  phantoms  composed  of  (a)  gelatin  and  (b)  De-ionized  water.  These  spectra  represent 
the  non-turbid  baselines  since  neither  phantom  contained  significant  scattering  media. 

Fig.  5.  LuTex  fluorescence  emission  experimentally  measured  at  different  source- 
detector  distances  in  homogeneous  scattering  gelatin  and  intralipid  based  phantoms.  The 


circular  domain  shown  represents  a  cross-section  of  the  cylindrical  phantom,  illustrating 
the  input  (red  arrow)  and  output  measurement  sites  (blue  arrows).  For  illustration 
purposes,  the  intensity  (shown  as  logarithm  of  intensity)  of  the  diffuse  excitation  field  is 
plotted  in  the  circular  region. 

Fig.  6.  Fluorescence  emission  spectra  through  a  liquid  phantom  of  DI  water  and  varying 
concentrations  of  intralipid  are  shown.  Only  spectra  measured  through  intralipid 
concentrations  of  0.1%  and  above  are  shown  for  the  two  detectors  farthest  from  the 
source. 

Fig.  7.  Absorption  coefficients  of  100%  water  and  300  nM  LuTex  are  plotted  for  discrete 
wavelengths  (green  line  with  data  points)  with  LuTex  fluorescence  emission  (red  dotted 
line)  in  dilute  solution,  indicating  the  overlap  in  absorption  and  emission  at  the  shorter 
wavelengths  from  700  to  750  nm. 

Fig.  8.  Diffusion  based  modeling  of  the  fluorescence  peak  through  a  60  mm  diameter 
phantom,  similar  to  Fig.  5.  The  dilute  emission  spectrum  of  LuTex  is  shown  as  dotted 
lines,  and  the  calculated  emission  spectrum  from  diffuse  emission  through  the  region  is 
shown  at  discrete  wavelengths  (black  circles).  Experimental  measurements  are  shown  in 
blue  for  intralipid  liquid  phantoms.  The  circular  domain  shown  represents  a  cross-section 
of  the  cylindrical  phantom  upon  which  the  intensity  (logarithm)  of  the  excitation  field 
calculated  from  the  diffusion  equation  is  plotted  for  illustration  purposes. 


Fig.  9.  Extinction  coefficients  of  chromophores  used  in  the  simulation  are  shown  in  (a). 
These  are  used  to  calculate  tissue  absorption  coefficients  at  any  wavelength  within  the 
NIR.  The  simulated  domain  contains  contrasts  in  the  various  absorbing  constituents,  the 
spatial  distributions  of  which  are  shown  in  (b). 

Fig.  10.  Fluorescence  yield  of  LuTex  is  calculated  from  the  drug  concentration  and  the 
fluorescence  quantum  yield.  Based  on  the  concentration  of  LuTex  assumed  in  Fig.  9,  the 
target  fluorescence  yield  is  shown  in  (a).  Reconstructing  images  from  data  collected  with 
long-pass  filtering  can  be  an  intractable  problem  as  shown  in  (b),  while  resolving  the 
emitted  spectrum  allows  accurate  recovery  of  the  true  fluorescence  activity,  as  shown  in 
(c).  All  images  shown  are  in  units  of  mm'1. 
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